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ABSTRACT 

ANALYSES OF QUASI-ISOTROPIC COMPOSITE PLATES UNDER QUASI-STATIC 
POINT LOADS SIMULATING LOW-VELOCITY IMPACT PHENOMENA 

Ajit Dhundiraj Kelkar 
Old Dominion University, 1985 
Director: Dr. Ram Prabhakaran 

Co-Director: Dr. W. Elber 

Composite laminates have high strength to density ratios that make 
them attractive for use in aircraft structures. However, the damage 
tolerance of these materials is limited because they have very low 
ultimate strains, no plastic deformation range, and no usable strength 
in the thickness direction. These limitations are very obvious when 
laminates are subjected to impact loads. Due to these impact loads, 
laminates suffer visible and invisible damage. To improve the material 
performance in impact requires a better understanding of the deformation 
and damage mechanics under impact type loads. 

In thin composite laminates, the first level of visible damage 
occurs on the back face and is called "back face spalling." A plate- 
membrane coupling model, and a finite element model to analyze the large 
deformation behavior of eight-ply quasi-isotropic circular composite 
plates under impact type point loads are developed. The back face 
spalling phenomenon in thin composite plates is explained by using the 
plate-membrane coupling model and the finite element model in 
conjunction with the fracture mechanics principles. The experimental 
results verifying these models are presented. The study resulted in the 



following conclusions: 


1. The large deformation behavior of circular isotropic membranes 
subjected to arbitrary axi symmetric loading can be obtained by 
solving a single nonlinear governing equation in terms of a 
radial stress. 

2. Accurate large deflection behavior of circular quasi-isotropic 
T300/5208 laminates can be obtained by using a simple plate- 
membrane coupling model. 

3. The functional form of deformed shape of the plate undergoing 
large deformations is different from the small deflection plate 
sol ution. 

4. The back face spalling action in thin composite laminates is a 
spontaneous action and can be predicted by using the fracture 
mechanics principles. 

5. Mixed mode (I and II) type deformations probably occur during 
back face spalling, however, mode I appears to govern the 
delamination growth during the spalling action. 
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Chapter 1 


INTRODUCTION 
1.1 Motivation 

Composite materials are being used increasingly in aircraft 
structures due to their high specific strength and stiffnesses, and the 
resultant weight savings. Traditionally, use of graphite fiber 

composites have been confined to secondary structures. However, 
requirements for reduced structural weight, improved aircraft 

performance, and efficiency are making composite materials increasingly 
competitive for expanded usage in primary, load carrying structures. In 
comparison to conventional metals, an understanding of the complex 

behavior of composites is still in its preliminary stage, and 
applications are based on knowledge gained through extensive 

experimental programs. Past experiences and experiments have confirmed 
that graphite fiber composite laminates have very low ultimate strains, 
no plastic deformation range, and no usable strength in the thickness 
direction. These limitations are very obvious when laminates are 

subjected to impact loads. Therefore the study of impact damage 

susceptibility of these composite laminates is increasingly important. 

Resin matrix composites are basically brittle materials, and the 
damage caused by impact differs from the damage on ductile metal 

structures. The ductile metals tend to develop indentations which are 
normally visible. On the other hand brittle materials like composites 
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tend to have both visible and invisible damage. Such damage is usually 
in the form of delaminations, matrix cracks and possible broken 
fibers. It can cause significant strength losses in composite materials 
[1-7].* The potential severity of visible and invisible damage has 

instigated numerous investigations into the impact behavior of 
composites. 

Past studies [8-12] indicate that in thin composite laminates first 
level of visible damage occurs on the back face of the laminate and is 
called as "back face spalling." A typical back face spalled laminate is 
shown in Fig. 1.1. Elber [8] conducted a series of tests on circular 
composite plates. He identified the sequence in which damage occurs in 
thin 8 ply T300/5203 quasi-isotropic graphite/epoxy circular plates 
subjected to center point loads. He showed that, first visible damage 
occurs in the lowest ply (8th ply) in the form of a matrix cracking 
parallel to the fibers. Further increase in the load levels results 
into the two dominant cracks in the eighth ply and delamina tins between 
plies 7 and 8 surrounds these two cracks. Additional increase in the 
load levels cause these delaminations and cracking to propagate unstably 
in the bottom ply in the fiber direction. This unstable growth of 
delaminations and propagation of the two dominant cracks result in the 
massive back face spalling action. During the back face spalling 
action, curved trajectory of the 8th ply center strip (formed by the two 
dominant cracks in the 8th ply) running under the load point in the 
fiber direction changes virtually to a straight line between the load 
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Fig. 1.1 Back Face Spalling in a Composite Plate 



4 


point and the support (Fig. 1.1). Several investigators [9-12] have 
also reported back face spalling damages in composite laminates under 
low velocity impact type point loads. V/hile considerable experimental 
observations have been made on the back face spalling, a theoretical 
understanding of the deformation mechanics and the back face spalling 
phenomenon is not fully known. The present " study was undertaken to 
develop analytical models, to understand the deformation mechanics, and 
to explain the back face spalling phenomenon, in thin composite 
laminates. These models will be particularly helpful to quantify the 
impact mechanics and to produce composites with improved resistance to 
the back face spalling and impact damage. 

1.2 Scope and Objectives 

Impact events are characterized into three velocity domains: 

1. High velocity or ballistic impact where the velocity, v > 600 
m/sec. 

2. Intermediate velocity impact (30 m/sec < v < 600 m/sec), such 
as foreign object damage (F.O.D.) impact on turbine fan blades, 
ground debris thrown up by wheels, etc. 

3. Low velocity impacts (v < 30 m/sec) are represented by handling 
damage such as dropping tools, foot steps and similar events. 

In this study, the effect of low velocity impact on composites is 
addressed. Many low velocity impact situations that an aircraft 
component is exposed to include: dropped tools, runway stones, tire 

blow out debris, ground collisions, etc. Due to low velocity impact. 
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composite laminates suffer both visible and invisible damage. To 
improve the material performance in impact- requires a better 
understanding of the deformation and damage mechanics under low velocity 
type impact loading. Therefore, the specific objectives of this study 
are: 

1. To develop an analytical model to understand the deformation 
mechanics in thin composite laminates under low velocity type 
impact loading. 

2. To develop an analytical model to expain the back face spalling 
phenomenon in thin composite laminates under low velocity type 
impact loading. 

3. To verify these analytical models by conducting experiments on 
thin composite laminates. 

For both the analytical and the experimental purpose quasi- 
isotropic laminates of T300/5208 graphite/epoxy material with stacking 
sequence of [45/0/-45/ 90] s are considered in this study. 

1.3 Review of Earlier Work 

The subject of impact related phenomena has been studied by many 
investigators utilizing many different approaches. Some of this work 
has been related to ballistic type impact, and hence, is not applicable 
here. The past studies indicate that considerable experimental 
investigations have been made to understand the low velocity impact 
problem in composite laminates. Literature review also indicates that 
theoretical investigations of the low velocity impact in composites are 
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1.3.1 Experimental Studies 

McQuillen and Gause [13] first studied the response of composite 
structural elements to low speed, transverse impact. They conducted 
series of experiments on graphite/epoxy composite laminates 76 mm wide, 
with a span 152 mm and with thickness of 2.13 mm. The laminates had 
(± 45/02/+45) stacking sequences and were simply supported on two 

edges. McQuillen and Gause concluded that inplane stress wave effects 
in the low velocity regime are negligible. They observed that the 
dynamic failure mechanism is identical to the static failure 
mechanism. Rhodes et al. [12] investigated the effect of low velocity 
impact damage on the compressive strength of graphite/epoxy panels. 
They concluded that the compressive strength of graphite-epoxy 
components can be reduced significantly by the effects of low velocity 
impact damage. Similar observations were made by Gause and Huang 

[14]. They observed significant reduction in the static strength of 
graphite/epoxy panels, when panels were subjected to low velocity impact 
loads. They examined the damage mechanisms in graphite/epoxy panels 
under dynamic and static loading. They concluded that in graphite/epoxy 
panels the static and dynamic damage mechanisms are identical. Card and 
Rhodes [15] studied the effects of low velocity impact on the compres- 
sive strength of graphite/epoxy structures. They conducted extensive 
tests on sandwich beams, laminated plates and stiffened panels. They 
studied the damage patterns and suggested possible sequence for laminate 
damage due to low velocity impact and found that, there was no visible 
damage on the impact surface, however, the side opposite to impact 
surface suffered a visible damage in the form of back face spalling. 
Their observations were similar to the earlier observations made by 
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Cristescu et al . [16] with the 0-90° ply fiberglass/epoxy composite 

plates. 

Gause [17] performed low speed, hard object impact on thick 
graphite/epoxy plates. He concluded that inplane wave effects are not 
important in low velocity impact on thick graphite/epoxy plates. He 
found that in thick composite plates the dominant damage mode is shear 
failure mode. He performed quasi-static tests on thick laminates and 
found that the damage mechanics in impact and quasi-static loading are 
identical. Gause et al. [18] studied the effect of low-velocity impact 
damage on the composte wing box. They identified that the first visible 
damage in composite skins occurs on the back surface (surface opposite 
to the impact point) in the form of spalling. Similar observations were 
made by Ramkumar [5], Bhatia [10], Hertzberg et al. [11]. 

1.3.2 Theoretical Studies 

As pointed out earlier only a few theoretical investigations were 
made. Some of the earlier work is discussed here. A theoretical model 
to study the damage and deformation mechanics due to low velocity impact 
problem in composites was developed by Horen [19]. He developed a 
theoretical model for a rectangular composite plate subjected to a 
quasi-static load equivalent to the impact load. Dobyns and Porter [20] 
and Dobyns [21] and Hayes and Rybicki [22] performed similar analysis 
for simply supported composite plates and Greszczuk [23] performed a 
quasi-static analysis for circular composite plates. However, all these 
analyses were based on the linear plate bending theory. 

Recent studies by Bostaph and Elber [24], Elber [8], Lai [25] and 
Shivakumar el al. [26], have shown that the laminates undergo large 
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deformation when the impact occurs at low velocity. Hence, it is 
necessary to develop- a more general analysis that includes large 
deflection effects. These types of analyses will help to predict 
accurately the deformations and failures in composite plates subjected 
to low velocity type impact loading. 

Bostaph and Elber [24] performed quasi-static impact tests on thin 
composite laminates. They used superposition approach to match the 
plate bending solution and a membrane solution under the load points. 
By using this superposition approach they obtained load-deflection 
behavior of circular composite plates undergoing large deformations. 
Shivakumar et al. [26] analyzed the impact damage problem using the Ritz 
technique and the large deformation theory for the circular thin quasi - 
isotropic laminates. The classical solution method used by Shivakumar 
et al. [26], Timoshenko [27], Washizu [28], to obtain the large 

deformation solutions of thin plates under the point loads is based on 
the Ritz technique, and only represents a membrane correction to the 

plate stiffness. Further the classical solution assumes that the 

functional form of the deformed shape of the plate is identical to the 
functional form of the deformed shape of the plate determined by the 

small deflection solution. The functional form of the deformed shape of 
the plate, however, is not identical for different load levels. Because 
at higher load levels the plate undergoes large deformations and the 
external load is partly equi 1 ibriated by the membrane action, the 
deformed shape of the plate can not be identical at all load levels. 

Accurate prediction of the deformed shapes are necessary to obtain 
accurate stresses in the plate. The plate analysis based on the 
classical solution, which assumes the functional form of the deformed 
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shape to be identical for all load levels, would incorrectly predict the 
stresses. One of the objectives, therefore, was to develop simple 
analytical models to predict the accurate large deformation shapes of 
thin circular laminates under quasi-static point loads. To obtain 
accurate large deformation shapes of these thin circular laminates, a 
finite difference and a finite element model is proposed in this 
thesis. As discussed earlier, several investigators have shown that the 
damages in thin composite laminates due to low velocity impact loads are 
similiar to the damages due to equivalent quasi-static loads, therefore 
in both the models proposed here, low velocity impact loads are replaced 
by equivalent quasi-static loads. These models are briefly discussed in 
the following section. 


1.4 Proposed Methods of Analyses 

In the present work an attempt was made to obtain simple analytical 
models for the large deflection analysis of thin circular quasi- 
isotropic laminates under the quasi-static point loads. The first model 
proposed here is a plate-membrane coupling model. This model does not 
have the constraining assumption that the deformed shape of the plate 
should be identical to the deformed shape of the plate determined by the 
small deflection plate solution. The plate-membrane coupling model 
proposed here assumes that the effect of flexural anisotropy of the 
axially quasi-isotropic laminates is small and tne plate behaves 
flexurally isotropic. 

1.4.1 Plate-Membrane Coupling Model 


A thin plate undergoing large deformations can be decomposed into 
two problems. First, a plate with shear and flexural stiffnesses but no 
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mid-plane extensional stiffness. Second, a plate with extensional 
stiffness, but no shear and flexural stiffnesses (membrane). Solutions 
of these two problems when coupled together will yield complete 
knowledge of displacements in quasi-isotropic laminates undergoing large 
deformations . 

Solution of the first problem is the well known small deflection 
plate analysis [27]. The second problem which is the analysis of 
circular membrane subjected to axisymmetric loading is more complex in 
nature, as it involves geometrical nonlinearity. The problem of 
circular membrane subjected to surface and edge loads has been studied 
by many investigators. Hencky [29], Dickey [30], and Shaw and Perrone 
[31] determined the deflections of a uniformly loaded membrane. 
Goldberg and Pifko [32] and Weinitschke [33,34] employed power series 
approaches to obtain the solutions of annular membranes. In addition to 
power series method Weinitschke [33,34] presented an integral equation 
approach to the solution of annular membranes subjected to surface and 
edge loads. Callegari and Reiss [35] studied the axisymmetric 
deformations of circular membranes subjected to uniform normal pressure 
by using the shooting method. The literature review indicates that 
limited numerical solutions are available for the nonlinear memorane 
problems. Furthermore these solutions are limited to the circular 
membranes subjected to uniform load. 

In the present work a more general formulation for the analysis of 
circular isotropic membranes subjected to arbitrary axisymmetric loading 
is presented. In this formulation a single nonlinear governing equation 
in terms of radial stress is used. The nonlinear governing differential 
equation is replaced by a set of nonlinear algebraic equations using 



11 


finite difference technique. These nonlinear algebraic equations are 
numerically solved to obtain stresses and displacements. 

By using the small deflection plate solution and the large 
deformation membrane soltuion in conjunction with the plate-membrane 
coupling model, the large deformation solution of the clamped circular 
quasi-isotropic plate subjected to quasi-static center point load is 
obtained. 

The plate-membrane coupling model uses a finite difference 
technique. A second analytical model is proposed in this work, to 
obtain the large deformation solutions of circular laminates uses a 
finite element method. The finite element formulation presented here is 
based on an appropriate linearization of strain displacement relations 
and uses an iterative method of solution. The formulation includes the 
effects of geometric nonlinearity. 

To study the geometric nonlinear problems many investigators [36- 
44], used the energy approach, in which the nonlinear strain displace- 
ment relation is linearized by equivalent linearization technique. 
However, earlier investigators have ignored the effects of inplane 
displacements in the formulation. This discrepancy was pointed out by 
Prathap and Vardan [45], Sarma and Vardan [46], Prathap and Bashyam 
[47], and Prathap [48]. However, so far no attempt has been made to 
include the inplane displacements in the finite element formulation of 
the circular plate subjected to a center point load and undergoing the 
large deformations. 


The reason for ignoring the effects of mid-plane displacements by 
earlier investigators may have been due to complexity of the nonlinear 
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finite element formulation. In the present finite element formulation 
the discrepancy -of ignoring inplane displacements is removed and a 
complete finite element formulation which includes inplane displacements 
is presented. The large deformation solutions of clamped circular 
plates subjected to central point loads obtained by using the plate- 
membrane coupling model, and the finite element model, are compared. 
Experimental verification of the plate-membrane coupling model and the 
finite element model is also presented. 

The second objective of this research is to develop a theoretical 
model to study the back face spalling phenomenon in thin composite 
laminates under low-velocity impact type point loads. Elber [8] 
conducted series of tests on thin eight ply T300/5208 quasi-isotropic 
graphite/epoxy laminates. He showed that the damage mode, in actual 
impact tests at velocities around 5 m/s, was very similar to the damage 
made in the static tests. He postulated that the massive back face 
spalling which often is the first visible sign of damage in laminates is 
related to the low peel-mode fracture toughness in brittle resins and 
that it is trackable by static analysis. In the present work, this 
phenomenon of the back face spalling is explained by developing a simple 
analytical model based on the quasi-static analysis. The model uses the 
large deformation shapes of the circular laminates in conjunction with 
the fracture mechanics principles. This model is verified by conducting 
the back face spalling experiments on thin composite laminates. 

1.5 Layout of Presentation 

In Chap. 1, the problem of low velocity impact in graphite/epoxy 
composites is discussed. The phenomenological aspects of the problem. 
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such as the large deformation behavior, the back face spalling etc. are 
introduced. The relevant literature on the study of the low velocity 
impact problem emphasizing the deformation mechanics and the back face 
spalling phenomenon, is reviewed. Against this background the present 
methods of analyses to study the deformation mechanics and the back face 
spalling phenomenon in thin composite laminates under quasi-static point 
loads are proposed. 

Chapter 2 presents a small deflection plate analysis for a clamped 
circular quasi-isotropic plate. A method for modeling flexurally 
anisotropic quasi-isotropic plate as an isotropic plate having the 
flexural stiffness components equivalent to the flexural stiffness 
components of the quasi-isotropic plate is presented. A finite 
difference method for computing displacements is described. The 
numerical results for the plate deformations, obtained by using the 
finite difference method are compared with the exact solution. 

In Chap. 3, a general analysis of circular isotropic membranes with 
clamped peripheral edges, subjected to arbitrary axiymmetric loading is 
presented. A single governing equation in terms of radial stress is 
derived. The solution of this nonlinear governing equation is presented 
by using the finite difference method in conjunction with Newton-Raphson 
method. 

Chapter 4 presents the plate-membrane coupling model to study the 
large deformation behavior of clamped circular plates. The model uses 
the plate solution (Chap. 2) and the membrane solution (Chap. 3), in 
conjunction with a coupling principle. The numerical results obtained 
by using the plate-membrane coupling model for clamped circular 
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laminates undergoing the large deformations are presented. These 
numerical results are compared with the classical solution. 

In Chap. 5, a finite element model which uses appropriate 
linearization of strain-displacement relations is presented. The model 
is based on a variational technique and uses the principle of minimum 
potential energy. The finite element model is then used to obtain the 
large deformation solution of clamped circular laminates subjected to 
central point loads. The solutions obtained by the finite element model 
and the plate-membrane coupling model are compared. 

In Chap. 6, plate mechanics experiments are described. The 
experimentally obtained large deformation shapes and the load-deflection 
curves of thin quasi-isotropic clampled circular laminates are compared 
with those obtained by using the plate-membrane coupling model (Chap. 4) 
and the finite element model (Chap. 5). 

Chapter 7 presents the back face spalling model. The back face 
spalling model is developed by using the large deformation shapes of 
thin circular laminates in conjunction with fracture mechanics 
principles. The analytical back face spalling model is verified by 
conducting back face spalling experiments on thin composite laminates. 
Several alternative approaches to reduce the back face spalling in 
composite plates are discussed. 

Finally, in Chap. 8, a summary of the highlights of the present 
work and some possible directions for further study of low velocity 
impact in composites are presented. 



Chapter 2 


SMALL DEFLECTION PLATE ANALYSIS 
2.1 Introduction 

In Chap. 1, a plate-membrane coupling model was introduced to 
obtain the large deformation behavior of quasi-isotropic laminates. The 
pi ate -membrane coupling model requires a small deflection plate solution 
and a large deflection membrane solution. In this chapter the small 
deflection plate analysis for a clamped circular quasi-isotropic plate 
is presented. To develop a simple one dimensional analytical model, 

quasi-isotropic plate is modeled as an isotropic plate with the flexural 
stiffness components equivalent to the flexural stiffness components of 
the quasi-isotropic plate. A finite difference method for computing 
displacements is described. A study to establish the convergence 
characteristics of the finite difference method is then presented. The 
numerical results for the plate deformations, obtained by using the 
finite difference method are compared with the exact solution. 

2.2 Plate Configuration 

Figure 2.1 shows an axisymmetrically loaded clamped circular plate, 
with thickness h and radius a. The plate was a quasi-isotropic 

laminate of T300/5208 graphite/epoxy material with stacking sequence 
[ 45/0/-45/90 J s . The material properties of the laminate are given in 

Table 2.1. 
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2.3 Assumptions 

The small deflection plate theory, generally attributed to 
Kirchhoff and Love [49], is used here with the following assumptions: 

1. The plate is initially flat, i.e. the plate has no initial 
curvature. 

2. The thickness of the plate is small compared to its other 
dimensions. 

3. The deflections are small compared to the plated thickness. 
(Maximum deflection-to-plate thickness ratios of 1/15 to 1/10 
are considered small) 

4. The slopes of the deflected middle surface are small compared 
to unity. (This follows from the assumption 3 above) 

5. The deformations are such that, straight lines initially normal 
to the middle surface, remain straight and normal to the middle 
surface, (i.e. transverse shear strains are not permitted) 

6. The material has the same elastic modulus in tension and in 
compression. 

7. The component of stress normal to the midsurface, a , is 
negl igible 

8. The strains in the middle surface, produced by inplane forces, 
are neglected. 

9. The plate is assumed to be flexural ly isotropic, even through a 
quasi-isotropic plate is anisotropic in flexure. 
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2.4 The Governing Differential Equation 

In this section the governing differential equation for axisym- 
metrical bending of a clamped circular plate subjected to arbitrary 
axisymmetric loading and undergoing small deformation is presented. 

Figure 2.2 shows an inf ini tisimal element abed cut out from the 
plate by two cylindrical sections ab and cd and by two radial sections 
ad and be. Consider the equilibrium of an element abed. The couple 
acting on the side cd of the element is 

M r r d0 (2.1) 

where M r is the radial moment per unit length. The corresponding couple 
on the side ab is 


(M r + 

r 

dr) (r 

+ dr) d0 

(2.2) 

Couples on the sides 

ad 

and 

be of the element are each 

M 0 dr 


where M 0 is the tangential moment per unit length. The components of 

d 8 

these couples in plane rz are each M Q dr sin (— ^) and for small 
d9, sin (-|.) = ^ . Therefore, these couples give a resultant couple 
in the plane rz equal to 


M 0 dr d0 (2.3) 

Denoting Q(r) as the shearing force per unit length of the 
cylindrical section of radius r, the total shearing force acting on the 
side cd of the element is Q (r) r d0, and the corresponding force on 
the side ab is [Q (r) + dr] (r + dr) d0, these forces give a 

couple in the rz plane equal to 
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Q(r) r d9 ^ + (Q (r) + (r + dr) d9 ^ (2.4) 

Neglecting the third and higher order terms in Eq. (2.4) the 
resultant couple in the rz plane equal to 


Q (r) r d0 dr 


(2.5) 


Summing up the moments, the equation of equilibrium of the element 

abed: 

dM 

(M + -g— dr) (r+dr) d9 - M r r d9 - M Q dr d9 + Q(r) r d9 dr = 0 (2.6) 


neglecting the third and higher order terms, one has, 

dM 

M r + W~ r ‘ M e + Q{r) r = 0 (2.7) 


Relations between curvatures and moments for pure bending, in 
cylindrical co-ordinate system are given by [27], 

£ 

dr' 


M = -D [i-£ + 1 %L] 


( 2 . 8 ) 


= - D [I 

0 


1 dw + v 


r <Tr 


(2.9) 


dr 


where D is a flexural modulus of the isotropic plate and is given by 

.3 


D = 


Eh* 

12 (1-v 2 ) 


( 2 . 10 ) 


where E is the Young's Modulus and v is the Poisson's ratio of the 
isotropic plate. 


Substituting Eqs. (2.8) and (2.9) into Eq. (2.7), the governing 
differential equation for the circular plate subjected to axisymmetric 
loading is 
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3 2 

d W + 1 d w 
dr^ r dr^ 


1 dw Q(r ) 
~P- "Hr D 


( 2 . 11 ) 


Equation (2.11) is a governing differential equation for the 
circular isotropic plate with the flexural modulus D and subjected to 
arbitrary axisymmetric loading. Denoting D eq as a flexural modulus of 
the quasi-isotropic circular plate the governing differential Eq. (2.11) 
for the quasi-isotropic circular plate was written as: 


where Q(r) 
obtained as 


3 2 

d w 1 d w 

1 ? 


1 dw Q(r ) 
°eq 


( 2 . 12 ) 


is a shear force per unit length at any radius r and can be 
(Fig. 2.3) 


2 nr Q(r) = 2 n l p (£) dC. (2.13) 

where p(£) is the intensity of loading at any radius 5. D eq in Eq. 
(2.12), is an equivalent flexural modulus for the quasi-isotropic 
laminate. The equivalent modulus D eq , can be obtained by equating 

bending energy of the quasi-isotropic laminate and an equivalent 
isotropic plate. A detail derivation to determine the equivalent 
flexural modulus D eq , is presented in the next section. 

2.5 Derivation of the Equivalent Flexural Modulus D eq 

Laminate extensional, coupling, and bending stiffnesses were 
calculated by using classical laminate plate theory [50] for a [45/0/- 
45/90 ] s quasi-isotropic laminate with elastic properties given in Table 
2.1. The extensional stiffness is independent of the polar angle 9 
(see Fig. 2.1) and is constant over the entire plate domain. This 
indicates that the laminate behaves perfectly isotropic for membrane and 
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inplane loadings. 

Since the laminate is symmetric, the coefficients of the coupling 
stiffness are found to be zero. The laminate bending stiffness depends 
on 0. However, an equivalent bending stiffness which is independent 
of 9 can be determined by equating the flexural strain energies of the 
clamped circular quasi-isotropic laminate and an equivalent circular 
isotropic plate. 

For a circular quasi-isotropic laminate the total strain energy in 
bending is 

U b = 7 J o U % {M}T {k} r dr d0 (2.14) 

where 

( M r 

(M> = jM e 

( f Ve 
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U K = J 2n /* {*c } T [D] {<} r dr d0 


(2.17) 


Equation (2.17) can be expanded as 


U b = 7 4* /o CK r D 11 + 2 °io + 2 D, 


r 0 w 12 


V 0 16 


+ 2 x p0 < e D 26 + Kg D 22 

+ k 2 q Dgg] r dr d9 (2.18) 

The coefficients D^, Dj^ , D 16* D 26» and D 66 in the 8c l- ( 2 -* 8 ) are 
all known functions of 6 [50]. Assuming that the curvatures are 

functions of r alone, then the integration over 0 in Eq. (2.18) can be 
carried out and Eq. (2.18) can be written as 

U b = ~7 4 ^ K r D 11 + 2 K r K e °12 + 2 *r0 K 0 D 16 


+ 2 < r g Kg D 2 g + Kg D 22 

+ 4q D 66 ] r dr (2 ' 19) 

where 

D ij = 7i l\' °ij d0; = 2 ’ and 6 

The bending energy of an isotropic plate with a flexural modulus 

Don. a Poisson's ratio of v and curvatures, k , k q , and k » (same 

eq rdrw 

as that of the quasi-isotropic laminates) can be written as 

U b (isotropic plate) = \ f* / 2u [k 2 D eq + 2 K p k q v 0q D gq 
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k 2 D 
6 eq 


+ K 


r0 


(l^H) D eq ] r dr de 


( 2 . 20 ) 


The flexural modulus D g q and the Poisson's ratio v are 
independent of 6. Hence integrations on 6 can be carried out and 
Eq. (2.20) can be written as. 


(isotropic plate) = — ^ 


2,1 / a U 2 D 


+ 2 < v D 
r eq r 0 eq eq 


1-v 


+ K e °eq + (-2 s5 ) V r dr (2 - 21 ’ 


For the same arbitrary curvatures, >c , k and k . on the composite 

r o ro 

laminate and the isotropic plate, the bending strain energies must be 
the same because of the assumed energy equivalence. Hence comparing 
curvature coefficient of Eqs. (2.19) and (2.21) equivalent bending 
stiffness D gq and v Poisson's ratio can be determined as. 


D 




de 


( 2 . 22 ) 


v 

eq 




d9 


(2.23) 






de 


(2.24) 


1-v 

(-tS2) 


eq 


1 

Tk 


t D 


66 


d6 


(2.25) 


If circular quasi-isotropic laminates were flexurally isotropic. 


then integrals 


26 


^ 0 16 de and ^ jf D 26 d9 (2.26) 

would be equal to zero. However, in the circular quasi-isotropic 
laminate under consideration, integrals represented by Eq. (2.26) were 
nonzero. However, for the material properties presented in Table 2.1, 
integrals were found to be about 2.5 percent of J° n d0 and 
hence neglected in this analysis. 

Equations (2.22) through (2.25) represent four equations, with two 
unknowns, D e q and v . The first two equations, Eqs. (2.22) and 
(2.23), are sufficient to evaluate the two unknowns D eq and . The 
third equation, Eq. (2.24), is essentially the same as Eq. (2.22) 
because of symmetries. For the material properties in Table 2.1, the 
value of D eq and v were found to be 5.688 N-m, and 0.31, 
respectively. These values of D e q and v eq , when substituted into 
the last equation, Eq. (2.25), satisfied the equation exactly. 

For those ply properties mentioned in Table 2.1, the value of 
equivalent Young's modulus E e q was obtained from an equivalent flexural 
modulus D e q as 







(2.27) 


and was found to be 53.3 GPa. This value of E e q is exactly the same as 
the inplane Young's modulus of the laminate obtained by using the ply 
properties in Table 2.1 and the classical laminate theory. 

2.6 Solution Method 


The solution of the clamped circular plate problem via the 



27 


classical method [27], is limited to relatively simple plate load 
conditions. In simple loading cases shear force Q (r) can be directly 
expressed as a function of r and since the governing equation, Eq. 2.12 
is Euler equation, it can be solved to obtain a closed form solution for 
the transverse displacement w. However, if loading conditions are more 
complex, the analysis becomes increasingly tedious. In such cases 
numerical methods are used to obtain the solution of the problem. Among 
the numerical techniques presently available, the finite difference 
method and the finite element method are most commonly used. In the 
present analysis the finite difference method is used. The solution of 
the plate problem using the finite element method is described later in 
Chap. 5. In the finite difference method, the derivatives in the 
governing differential equation are replaced by difference quantities at 
selected points, called nodes. The details of the finite difference 
method are as follows: 

The basic finite difference expressions follow logically from the 
fundamental rules of calculus [51]. For a continuous function w = f(r) 
(Fig. 2.4), the first, second and third derivatives of the displacement 
at any node n can be written by using the central difference method 
[52] as: 


dw _ w(n+l) - w(n-l) 
cF ln; * 2 ~^~aF 


(2.28) 


d w , . w(n+l) - 2 w(n) + w(n-l) 

— (n) = 2 

dr 2 (A r) 2 


(2.29) 


d 2 w , . w(n+2) - 2 w(n+l) + 2 w(n-l) - w(n-2) 

— (n) = 

dr 2 • (A r) J 


(2.30) 



ArArArArArAr 


Finite Differences 


for a Continuous Function 
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The governing differential equation is then transformed to an 
algebraic equation, by substituting the appropriate finite difference 
expressions. The solution of the differential equation thus reduces to 
the simultaneous solution of a set of linear, algebraic equations, 
written for every nodal point. 

2.6.1 Application of the Finite Difference Method for the Solution of 
Arbitrarily Loaded Clamped Circular Plate 

To solve the governing Eq. (2.12), the solution domain was 

discretized into m regions and (m + 1) nodes. Denoting w(n), (n), 

2 3 ar 

d w d w 

— 7 (n), and — -• (n) as the first, second, and third derivatives of 
dr dr J 

the displacement at the n node, the governing Eq. (2.12) was written as, 
at any node n: 


2 

d w 


(n) + 


? 

d w 


dr 


n dr 


(n) - (n) 



n pU) l dl 

r D 
n eq 


(2.31) 


First, second, and third derivatives of the transverse displacement 
w from the governing Eq. (2.31), were replaced by finite difference 
quantities given by Eqs. (2.28), (2.29), and (2.30). 

To simplify the evaluation of the integral in the above equation 
the following assumption was made. Consider an i region with 
r*i_l < r < r.. Although the applied load varies within the region 
r._j < r < r^, the load will be assumed to be uniform in this region 
with a value p^. The magnitude of p^ is assumed to be equal to the 
value of the load at the mid-point of this region, i.e., at r = (r^ + 
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r j _ 2 ) / 2. As the number of regions in the model become large, the size 
of each region decreases and hence the variation of the load within each 
region becomes insignificant. With this assumption the governing Eq. 
(2.31), reduces to 


pf‘"> + ^$ <n, ;r&< n) = 

n 


n 2 2 

I p (r - r ) 
i=l i i i-1 

2 r n D on 

n eq 


(2.32) 


2.6.2 Boundary Conditions 

For the circular plate with axi symmetric loading the boundary 
conditions were 

(a) Both the transverse displacement w and the slope dw/dr equal zero 
at the clamped edge (r = a). 

(b) The slope dw/dr equal zero at the center (r =0). 

Using the governing Eq. (2.32) in the form of finite difference 
quantities at each node and transforming the above boundary conditions 
into finite difference quantities, a set of simultaneous algebraic 
equations were obtained. This set contained number of unknowns that are 
equal to number of nodes in the solution domain. These algebraic 

equations were solved simultaneously, to obtain the transverse 
displacement at each node. To illustrate the method, a very simple 4 
region idealization was chosen. The details of the method are presented 
in the followng section, for this simple idealization. 

2.7 Illustrative Example 


Consider a clamped circular quasi -isotropic plate with radius a 
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and thickness h. The plate is subjected to a uniform loading over the 
region 0 < — < 0.5. The objective is to obtain the plate 

a 

displacements using the finite difference method. 

2.7.1 Problem Formulation 

The solution domain 0 < — < 1.0 is discretized into 4 regions and 

cl 

5 nodes (see Fig. 2.5). At the node 0, the boundary condition is slope 
dw 

= 0. This conditions can be transformed into the finite difference 
form as follows: 

By using Taylor's series expansion, the displacements at the nodes 
(1) and (2) can be written as 


w(l) = w(0) +^1^(0) + (0) (2.33) 

dr 

w(2) = w(0) + ij£ (0) + i2 f { )Z (0) (2.34) 

l. or z. dr * 

2 

d w 

From the Eqs. (2.33) and (2.34) by eliminating terms — (0) and using 

d d r ^ dw 

boundary condition, the slope equal zero at r = 0 or (0) = 0, 

Eq. (2.34) is written as: 

4 w(l) - w( 2) - 3 w( 0) = 0 ’ (2.35) 


The governing differential equation at any node n for the quasi- 
isotropic plate is (Eq. (2.32)) 


^5 (n) + -r {n) ' 4 W (n) = " 

r n? 7 W i= 1 

n 


P, (rf - rf_ x > 

2 r„ D 
n eq 
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Fig. 2.5 Four Region Idealization for the Plate Uniformly 
Loaded Over the Region 0 £ r/a < 0.5 
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For the node no. 1 the first and second derivatives can be 
expressed in the finite difference form, by using central difference 
Eqs. (2.28) and (2.29) as: 


dw i ■} \ _ w(2) - w(0) 

w (1} r~ r ~Ar 


(2.36) 


d w , lX _ w( 2) - 2 w(l) + w(0) 
— 7 l i) ? 

dr (A rr 


(2.37) 


.3 

Q W 

The third derivative i.e. — (1) can be obtained as follows: 

dr 

By using Taylor's series expansion the displacements at the nodes 
(2) and (3) are written as 


w<2) -w(l) + (1) (1) * dfw ( 1 ) 

i. or e.. dr <i j. 


(2.38) 


? ? 3 3 

n\ / . 2Ar dw / 1 \ . (2Ar) d w , (2Ar) d w M v ,, , n , 

w( 3) = w(l) + — jy- (1) + -■ 2 - | ^ ' 1 J + — jt 3 (D (2.39) 

dr ‘ dr 


Eliminating (1) from the Eqs. (2.38) and (2.39) the following 

equation is obtained 


2 3 

2 w( 2) - w(3) = w(l) - (Ar ) 2 1" (1) - (Ar ) 3 L* (l) 


(2.40) 


dr 


dr‘ 


But from the Eq. (2.37) 


d w ,,, w( 2) - 2 w(l) + w(0) 

— k- (1) - - O 

dr 2 (Ar / 
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Substituting Eq. (2.37) into Eq. (2.40), equation for 
obtained as follows: 



is 


- 3 w(2) + w(3) + 3 w(l) - w(0) = d J w 
(Ar)^ dr^ 


(1) 


(2.41) 


Using Eqs. (2.36), (2.37), and (2.41) in Eq. (2.32), at node no. 1, the 
governing Eq. (2.32) is 


0 ^2 r 2 (Ar) r 2 ( Ar ) 2 ( Ar ) ^ ^ 1 ^ r 1 (A r) 2 

+ w 2 [- ^ + 1 —2 ~ — 2-j] + w [— = 

2 2 r 2 (Ar) r x (Ar) 2 (Ar) J J (Ar) J 


(Ar) 


P < r l - r o ) 

2 r. D 
1 eq 


(2.42) 


Using the central differences (Eqs. (2.28)-{2.30) ) , at the node no. 
2, the governing Eq. (2.32) is. 


w [- j] + w [ ^ + 

2 (Ar) ° - ,A 


7 + t] + w ? [“ 

2 r^ (Ar) r 2 (Ar) 2 (A rr 


+ w 3 t* 


2 r^ (Ar) r 2 (A r)"^ (Ar)^ 


M + w 4 t 


1 


2 (Ar)’ 


r 2 (Ar)‘ 


•] = 


P tr| - r*) 

2 r~ D 
2 eq 


(2.43) 


and at the node no. 3, using Eqs. (2.28)— (2. 30 ) , the governing Eq. 
(2.32) is. 


W, [- i-J + w 2 [ ♦ * * -i-j] * W [ ] 

1 2 (ArP 2 (r 3 ) 2 Ar r 3 (Ar) 2 (ArP r 3 (Ar) 2 



35 


+ w 4 T + " — 2 " — ^-3] + w 5 t ^-3] 

4 2 r 3 (Ar) r 3 (Arp (ArP 5 2(ArP 


P ( r 2 ~ r 2 Q ) 


2 r, D 

3 eq 


(2.44) 


W 5 in the Eq. (2.44) is a displacement at fictitious node 5 (Fig. 

dw 

2.5). The boundary conditions at r = a are, (^-) = 0 and w = 0 and 
can be expressed in the finite difference form as 


dw /yn _ w( 5) - w( 3) _ n 

w rwi U 


(2.45) 


or 


w( 5) = w(3) (2.46) 


By using the Eq. (2.46), Eq. (2.44) is written as 


W 1 + w ? f + " — 7 + — ^ 3 ] 

2 ( Ar) J 2 r 3 (Ar) r 3 (Ar) (Ar) 

+ W 3 ~ — 2 + " — 3-1 + w 4 t T + " — 

r 3 (Arp 2 (Ar) J 2 r 3 (Ar) r 3 (Ar) 

P (r 2 - r 2 ) 

2 r, D Qr , 

3 eq 

Lastly at the node no. 4 the displacement w = 0. Therefore 

w( 4) = 0 (2.48) 


±-1) 


(Ar) 


(2.47) 


Thus Eqs. (2.35), (2.42), (2.43), (2.47), and (2.48) contains five 
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unknowns, w(0), w(l), w (2), w{3), and w(4) and can be solved 
simultaneously to obtain the unknown displacements. A simple example 
problem is presented here to illustrate the above method. 

Consider a clamped circluar plate of radius 25.4 mm and thickness 
1.05 mm. Let flexural modulus D e q, of the plate is 5.688 N.m. The 

displacements in the plate under uniform loading of intensity 100,000 
N/m 2 over the region 0 < r < 12.7 mm can be obtained as follows: 


Since the plate radius is 25.4 mm Ar is (25.4/4) = 6.35 mm. Also 
r l* r 2» r 3» anc * r 4 are 6.35 mm » 12.70 mm, 19.05 mm, and 25.4 mm 
respectively. Substitution of these values in Eqs. (2.35), (2.42), 

(2.43), (2.47), and (2.48), gives five simultaneous equations, which can 
be expressed in the matrix form as 


-3 

4 

1 

0 

0 


w(0) 1 

A 

0 

1 

2 

-5 

2 

0 


w( 1) / 

p ( Ar ; 

1 

-4 

13 

-8 

-5 

4 

■ 

w( 2) \ 

l) 

eq 

8 

0 

-9 

25 

-3 

-13 


w( 3) ( 

12 

0 

0 

0 

0 

1 


w(4) | 


0 


or in concise form 

[A] (w) = [F] 

where 


p = 100,000 N/m 2 


(2.49) 


(2.50) 


Ar = 6.35 mm, and D 


eq 


5.688 N-m 


From equation (2.50) (w) is 
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{w} = [A]" 1 [F] (2.51) 

Solving equation (2.51) for {w}, the nodal displacement are found to 
be: w(0) = 0.076226 mm, w(l) = 0.06677 mm, w(2) = 0.038410 mm, w(3) = 

0.005434 mm, and w(4) = 0.0 mm. 

2.8 Convergence Study 

To study the convergence of this method, the circular plate was 
idealized into m number of regions with (m + 1) nodes. A systematic 
convergence study was then made by doubling the number of regions. The 

number of regions m, used in this convergence study were 4, 8, 16, 32, 

60, and 64. Figures 2.6 and 2.7 present the relative errors in the 
maximum deflections normalized with respect to the exact central 
displacement, for the two different cases: case (i) uniformly loaded 

plate and case (ii) a plate loaded uniformly over the regions 
0 < j < 0.5. The analysis indicated that about 60 regions were 

p 

necessary for a plate loaded uniformly over the region 0 < — < 0.5 and 

a 

for a uniformly loaded plate to yield a solution which is within 0.1 
percent of the exact solution. 

2.9 Results and Discussion 

In this section the central displacements obtained by using the 

finite difference method are compared with the exact solution. Next 
deflected shapes of the plate obtained by using the finite difference 
method, are compared with the exact solution. 

The comparison of the central displacements obtained by using the 
finite difference method and the exact solution is shown in Table 2.2. 
The percentage errors in the central displacements obtained by using the 





Relative Error 


Fig . 2 



.7 Convergence Study for the Plate Uniformly Loaded 
Over the Region 0 _<_ r/a <_ 0 . 5 
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finite difference method relative to the exact solution were found to be 
0.055 percent and 0.072 percent for the uniformly loaded plate and for a 
plate loaded uniformly over the region 0 < — < 0.5 respectively. 

a 

Figures 2.8 and 2.9 show the percentage relative errors in the deflected 
shapes obtained by using the finite difference method and the exact 
solution, for the uniformly loaded plate and the plate loaded uniformly 
over the region 0 < — < 0.5. In both the cases maximum percentage 

a 

error is within 0.1 percent of the exact solution and occur at the 
center of the plate. 


2.10 Concluding Remarks 

In this chapter a numerical solution method to obtain the small 
deflection plate solution for a clamped circular quasi-isotropic plate 
was, developed. The quasi-isotropic plate was modeled as an isotropic 
plate having the flexural stiffness components equivalent to those of a 
quasi-isotropic plate. By using a finite difference method, the plate 

governing differential equation was replaced by a set of algebraic 
equations. These algebraic equations were solved simultaneously to 

obtain the transverse displacements for arbitrarily axisymmetrically 
loaded clamped circular quasi-isotropic laminates. The numerical 
solution was found to be within 0.1 percent of the exact solution for a 
uniformly loaded plate and for a plate loaded uniformly over the region 
0 < — < 0.5 by using 60 region idealization. 

cl 

In order to study the large deformation behavior of the quasi- 
isotropic plate, membrane effects should be incorporated in the small 
deformation soltuion. The following chapter presents the large 

deflection membrane analysis. 



Fig. 2.8 Comparison of the 
Using the Finite 
Deflected Shape f 
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Fig. 2.9 Comparison of the Deflected Shape Obtained by 

Using the Finite Difference Method and the Exact 
Deflected Shape for the Plate Uniformly Loaded 
Over the Region 0 < r/a < 0.5 


I 
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Table 2.1 - Elastic Properties of the Plate 


Modulus GPa 

Material E-q E 22 Poisson ratio 

Gr/Ep lamina 131.0 13.0 6.4 0.34 

Number of plies = 8 
Laminate thickness = 1.05 mm 
Stacking sequence: [45/0/-/45/90] s 

(Subscripts 11 and 22 correspond to the longitudinal and transverse 
direction of fiber) 
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Table 2.2 - Comparison of Maximum Center Displacements for 

Uniformly Loaded Plate and for a Plate Uniformly 
Loaded Over the Region 0 < — < 0.5 

cl 



Central transverse 
at r 

displacement w, 
= 0. 

Type of 
loading 

Present* 
sol uti on 

Exact 
sol ution 
(ref. 27) 

Uniform loading 

4 

0.015633 

eq 

4 

0.015625 

eq 

Uniform load over the 
region 0 < — < 0.5 

ct 

4 

0.009994 

eq 

4 

0.009987 

eq 

*0btained by using 60 regions 

idealization 




Chapter 3 


LARGE DEFLECTION MEMBRANE ANALYSIS 
3.1 Introduction 

In Chap. 1, a plate-membrane coupling model was proposed to obtain 
the large deformation behavior of quasi-isotropic laminates under point 
loads simulating low-velocity impact. This model requires a small 
deflection plate solution and a large deflection membrane solution. In 
the previous chapter a numerical method to obtain a small deformation 
solution for a clamped circular quasi-isotropic plate, subjected to 
arbitrary axisymmetric loading was developed. In this chapter a 
numerical method is developed to obtain a large deflection membrane 
solution. Again a quasi -isotropic circular laminate is considered in 
the membrane analysis. This quasi-isotropic laminate has only 
extensional stiffness but no shear and flexural stiffnesses. As pointed 
out in the previous chapter, quasi-isotropic T300/5208 circular 
laminates with the stacking sequence of [45/0/-45/90] $ are axially 
isotropic, i.e., extensional stiffness of these laminates is independent 
of the polar angle (Fig. 3.1). Hence, these circular laminates, with 
only extensional stiffness, can be modeled as circular isotropic 
membranes., In this chapter, a general analysis of circular isotropic 
membranes with clamped peripheral edges, subjected to arbitrary 
axisymmetric loading is presented. A single governing equation in terms 
of radial stress is developed. The solution of this nonlinear governing 
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Fig. 3.1 Membrane Configuration 
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equation is presented by using the finite difference method in 

ft 

conjunction with Newton-Raphson method. Three loading cases namely (a) 
uniformly loaded membrane, (b) a membrane with uniform load over an 
inner portion, and ( c ) a membrane with ring load are analyzed and the 
results are compared with the classical solution. 


3.2 Membrane Configuration and Strain-Displacement Relations 

Figure 3.1(a) shows an axi symmetrical ly loaded clamped membrane 
with thickness h and radius a. Since the deflection surface is 
axi symmetrical , the displacement can be resolved into two components: 
(1) a component u in the radial direction, and (2) a component w 
perpendicular to the plane of the membrane. From large deflection 
theory [27], strain-displacement relations for an isotropic membrane are 
obtained as follows: 

Consider an element AB of length dr as shown in Fig. 3.2. The 
radial strain e r from Fig. 3.2 can be expressed as 


(u 


£ = 


+ W * r) - 
dr 


ds - dr 
dr 


but 


(3.1) 


ds - dr , , 7 /dw»2 

3r - / 1 + ’dr’ 


1 



( dw,2 


(3.2) 


Substituting Eq. (3.2) into Eq. (3.1) 




( dw ) 2 
"dr 


(3.3) 


The tangential strain e 0 from the Fig. 3.3 is 
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^ _ (r + u) d6 - r d0 _ u 

6 RTS r 


(3.4) 


By using Eqs. (3.3) and (3.4), the radial and tangential stresses 


are: 


- _ L eq rdu 1 /-dw^ %1 u-i 

r " 7i“"2; W 7 W eq r - 1 
1 -V eq 


(3.5) 


a. = ^ fid + v (*1) 2 ] 


' e TIT 7 ) [? + ^ 77 

eq 


(3.6) 


where E e q and v are the equivalent Young's modulus and Poisson's 

ratio of the isotropic membrane and were calculated by using the 
technique presented in Chap. 2. 

Note that the assumption of large deflection but with small strains 
is made in this formulation. 


3.2.1 Equilibrium Equations 

The equation of equilibrium in the radial direction is obtained as 
follows: 

Consider stresses on an infinitesimal element abed of unit 
thickness as shown in Fig. 3.3. Summing the forces in the radial 
direction along the radial line (6 + one has 


da r 

( <?r + -g-p— dr) (r + dr) d0 - a p r d0 


dr sin 


d0 

~2 


(3.7) 


Since 

neglecting 


d0 is small, sin (^) may be replaced by (-^) . By 
higher order terms the equation of equilibrium in radial 
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direction is 


da r 

a r ~ a e + r W~ = 0 


(3.8) 


The equation of equilibrium in the direction perpendicular to the 
plane of the membrane, with general axisymmetric loading p(r) is 
obtained as follows: 

Consider a stretched circular membrane element with general 
axisymmetric loading p(r) as shown in the Fig. 3.4. Consider the free 
body equilibrium of the deflected membrane, i.e. the vertical component 
of the reaction must equal the total load. Therefore from Fig. 3.4 

2 n r h cr^ sin a - /£ p(£) 2 n £ d£ = 0 (3.9) 
dw 

For small a, sin a = tan a = 

Therefore the equation of equilibrium in the direction per- 
pendicular to the plane of the membrtane with general axisymmetric 
loading p(r) is 

2 n r h a r j! - p (^ 2 n ? = ° (3.10) 

The stress and strain displacement relations (Eqs. (3.3)-(3.6) ) and 
equilibrium Eq. ((3.8) and (3.10)) when combined form four nonlinear 
partial differential equations with four unknowns a , a Q , u, and w. 

To obtain a numerical solution of these governing equations the 
conventional approach is to reduce these four nonlinear equations to two 
equations containing only u and w displacements. This conventional 
approach is applicable when u and w displacements are of the same order 
of magnitude, for example, a case of uniformly loaded membrane as shown 
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by Shaw and Perrone [31], and Kao and Perrone [53]. However, the 
numerical solution of arbitrary axi symmetrically loaded membranes cannot 
be obtained by using the conventional approach because of the large 
differences in magnitudes of u and w displacements. Due to these large 
differences the two governing equations involving u and w displacements 
are difficult to satisfy. Hence an alternate formulation of the 
membrane problem is presented below. 


3.2.2 Derivation of a Governing Equation 

Using the stress and strain displacement relations (Eqs. (3.3)— 
(3.6)), the radial displacement u was expressed as 


u ' r- (o 8 ' V, a r> 

eq M 

du 


(3.11) 


The term was obtained by differentiating the radial displace- 

ment u, with respect to r. 


du 

dr 


i _ da R do 

- — ( a n - v a ) + 1-3 — - v —r-j 

. 0 . eq r t ^dr eq dr ; 

eq eq 


(3.12) 


The strain in the radial direction e r was expressed in terms of 
the radial and tangential stresses and by using the strain-displacement 
relations (Eqs. (3.3) and (3.4)), 


e = T — ( o - v oj = [■ 
r E r eq 0 L 
eq 


du , 1 fdw^-i 

W 7 W J J 


(3.13) 


du 


By using from Eq. (3.12) in equation (3.13) the following 

relationship was obtained:' 


{0 0 ' 


°r } < 


1 + 
~T 


eq- 


eq 


"eq 


r d0 0 

~ v eq 


d o 

r l 

^r~ J 


1 r dw > < 

7 [ W } 


(3.14) 


From the equation of equilibrium in radial direction (Eq. (3.8)) 
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u r 

(^p-) was written as 


d °r . 1 , 
dr r a r 


(3.15) 


Substituting Eq. (3.15) in Eq. (3.14) and rearranging the terms the 
following relationship was obtained: 


r , da e d V 1 /-dw^ _ Q 
{ W~ + Hr~' 7 [ W 1 0 


(3.16) 


Further, was obtained by differentiating the equation of 

equilibrium (Eq. (3.8)) in radial direction with respect to r, as 


do 9 , do r d2 ° r 


(3.17) 


— y 

Substituting this value of (-^p— ) in governing Eq. (3.16), the 

following relationship was obtained: 


r r _ u r . r ^ 1 r dw^2 n 

- (3 3r* r rr)*7y = 0 


(3.18) 


By using equation of equilibrium (Eq. (3.10)) in the direction 

( dw i 2 • 

perpendicular to the plane of membrane, (.pjp-j was written as 


r /q pU> 2 u i d^ i 2 

5 tt r- Ft ~rx J 


(3.19) 


Substituting this value of (^p)^ in the governing equation, Eq. 


(3.18) gives 
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d \ , 1 A PK> ? «, 2 . „ 
w T ' ur ) i *■ — nrr — J ~ 0 


- do 
S- (3 -zJL + r 


(3.20) 


"eq dr " r 

Equation (3.20) is a nonlinear differential equation in terms of 
the radial stress that governs the large deflection response of the 
membrane with arbitrary axisymmetric loading. Several investigators 
obtained the governing equation in terms of the radial stress a , 
similar to Eq. (3.20). Dickey [30] and Weinitschke [33] presented the 
differential equation for uniformly distributed load. Callegari and 
Reiss [35] obtained the differential equation for a membrane with an 
arbitrary axisymmetric loading. 


3.3 Solution Method 

To solve the governing equation (3.20) for the radial stress, 

Dickey [30] used integral equation method, Weinitschke [33] used 

integral equation and power series approaches and Callegari and Reiss 
[35] used the shooting method. In contrast, here a numerical method of 
solution is proposed. The nonlinear differential equation was replaced 
by a set of nonlinear algebraic equations using difference quotiets. 
Then using the Newton-Raphson method [52], these nonlinear algebraic 
equations were solved numerically to obtain the stresses and 

displacements. The details of the solution method are as follows: 

The solution domain was discretized into m regions and (m + 1) 
nodes. Denoting o r (n) as the radial stress at the nth node, the 

governing equation, Eq. (3.20) was rewritten as, at any node n: 



[3 


da r 

W~ 


(n) 


+ r 

n 




0 1 pU) 5 dS 2 

r h a (n) ^ ~ 0 

n n 


(3.21) 
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To simplify the evaluation of the integral in the above equation, 
the same assumption as discussed in previous chapter will be made. 
Consider an ith region with r._j, < r < r^. Although the applied load 

varies within the region r^ < r < r. , the load will be assumed to be 

uniform in this region with a value of p^ . The magnitude of p^ is 

assumed equal to the value of the load at the midpoint of this region, 

i.e., at r = (r^ + r^_^)/2. As the number of regions in the model 

become large, the size of each region reduces and hence the variation of 
the load within each region also becomes insignificant. With this 
assumption the governing equation, Eq. (3.21), reduces to 



+ 


1 

F 


da 

£ 3 IT <nl + r n 


dV 


dr‘ 


(n) ] 



1 

h ajb) 


n 

Z 

i=l 


Pj <r‘ - r 


2 

i-1 



(3.22) 


First and second derivatives of the radial stress from the 
governing Eq. (3.22) were replaced by using the central finite 
differences as: 


and 


da a (n+1) - a (n-1) 

dr - (n) = FTAr] 



(n) 


o. (n+1) - 2 a r (n) + a^ (n-1) 
(Ar) 2 



57 


3.3.1 Boundary Conditions 

For the circular isotropic membrane with axisymmetric loading the 
boundary conditions were: 

1. Both the radial displacement u and the transverse displacement 
w equal zero at the fixed edge (r = a). 

2. The radial displacement u and the slope equal zero at the 

center (r = 0). 

Since the governing equation, Eq. (3.20) was derived in terms of 
the radial stress, the transformation of boundary conditions was done by 
using stress and strain-displacement relations (Eqs. (3.3)-(3.6)) and 
equilibrium euquations, Eqs. (3.8) and (3.10). When u = 0 and r = a are 
substituted into equation (3.11) for radial displacement, the boundary 
condition no. 1 above can be transformed to 

°9 = v eq °r* at r = a (3.23) 

Using this boundary condition in the equilibrium equation in the radial 
direction (Eq. (3.8)), Eq. (3.23), was further transformed to 

da r 

°r {1_v eq ) + r tfr" =0 ’ atr = a (3.24) 

The boundary condition (2) above is a statement of symmetry about r = 0, 
the center of the membrane. This symmetry condition can be expressed in 
terms of the radial stresses, 

do r 

=0. at r = 0 (3.25) 

Using the governing equation, Eq. (3.22), 
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[3 


eq 


da r 

W~ 


(n) + r 


d 2 o. 


dr 


(n) ] 


1 

S 


n 

{ 2 

i=l 






n 


° r ( n ) 


} = 0 


and boundary conditions (Eqs. (3.24) and (3.25)) in the form of finite 
difference quotients at each node, m + 1, nonlinear algebraic equations 
were obtained. These algebraic equations contained m + 1 unknowns, 
viz., o r (0) ... a r (m) . By using the Newton-Raphson technique as 
described in Appendix A, these (m+1) equations were solved to obtain the 
radial stress at each node. 

do r 

Once the radial stresses were known, at each node, the -s — at 

dr 

each node was calculated by using the finite difference method. Using 

da r 

these values of a r and in Eq. (3.8), the tangential stress a Q 

at each node was obtained. By using Eq. (3.11), the radial displacement 

dw 

u at each node, was obtained. The slope pp at any node was obtained 
using Eq. (3.19). Then using the boundary conditions. 


w = 0 at r = a. 


^ = 0 at r = 0, 

dw 

and the finite difference representation of the slope, pp , at each 
node, (m + 1) simultaneous equations with w(0) . . . .w(m) unknowns, 
were obtained. These (m + 1) equations were then solved to obtain 
transverse displacement w at each node. The details for a simple four 
region idealization to obtain radial stresses at each node are presented 
next: 

3.3.2 Illustrative Example 


Consider a clamped circular membrane with radius a and thickness h. 
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The membrane is subjected to a uniform loading over the region 
0 < — < 0.5. The objective is to obtain the membrane displacements 

a 

using the finite difference method. 

The solution domain 0 < £■ < 1.0 is discretized into 4 regions and 

do r 

5 nodes (Fig. 3.5). At the node 0, the boundary condition is = 0 
(Eq. (3.25)). This boundary condition can be written in the finite 
difference form using Eq. (2.35) as 


4 a r (l) - 3 a r ( 0) - a^l) = 0 


(3.26) 


The governing differential equation at any node n for the membrane 
is (Eq. (3.22)) 


d 2 «. 


_ - d 2 o . n p . (r? - r? ) 2 

_ [3 („) r„ -j- <») ] j U ■ r 1 *0 


For the nodes 1, 2 and 3 the first and second derivatives can be 
expressed in the finite difference form, and by using central difference 
equations, Eqs. (2.28) and (2.29), Eq. (3.22) is 
At the node 1: 

- 0.5 a.(0) o 2 (l) - 2 oJ( 1) + 2.5 0.(2) o 2 (l) 


p 2 a 2 E 


eq 


128 h 


(3.27) 


At the node 2: 


o.(D 


( 2 ) 


8 o^(2) 


+ 7 


° r (3) 


o 2 (2) 

r 


2 2 F 
p a E 


eq 


32 h 


(3.28) 
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At the node 3: 


1.5 a (2) of (3) 
r r 


6 a^(3) + 4.5 o p (4) a 2 (3) 


(3.29) 


2 2 
P E € 

216 h 1 


At the node 4, the boundary condition (Eq. 3.24)) needs to be 
satisfied. This condition is 


o r <4) (l-» > * r 4 -gJ: (4) = 0 


(3.30) 


The -j— - (4) in the Eq. (3.30) can be obtained as follows: 

Using the Taylor's series expansion c r (3) and 0,(2) can be 


written as: 


2 

z A , do i . 1 2 do 

(3) = o (4) (4) +i£l ^ (4) 

dr 


(3.31) 


<Y (2) = °r (4 > -TT^ 1 (4) + ^TT “T" (4) 


(3.32) 


By eliminating — (4) terms from the Eqs. (3.31) and (3.32) and 
do r dr^ do. 

solving for (4), the (4) is obtained as follows: 


4 o r (3) + o^(2) + 3 o r (4) da^ 


(3.33) 


Substituting Eq. (3.33) into Eq. (3.30) the equation at the node 4 is 
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2 0 r (2) - 8 0 p (3) + a,(4) (6 + (1 - v eq )) = 0 (3.34) 

Equation (3.26), (3.27), (3.28), (3.29), and (3.34) are the five 

nonlinear algebraic equations and they contain five unknowns o r (0), 

0 r (l)» 0 p (2), 0 r (3), and 0 r (4). These nonlinear algebraic equations 

are solved by using Newton-Raphson' s method (see Appendix A) to obtain 

the unknown radial stress at each node. As discussed earlier, once the 

radial stresses were known at each node, the equation of equilibrium in 

the radial direction (Eq. (3.8)), was used to obtain the tangential 

stress, o Q , at each node. The radial displacement, u, and the slope, 
rl w 

, at any node were obtained by using Eqs. (3.11) and (3.18), 
respectively, and the o r and o Q values at that node. Then using the 
boundary conditions. 


w = 0 at r = a. 



at r 


0 


and the finite difference representation of the slope, -g^- , at all the 
nodes, the transverse displacements, w, were determined. 


3.4 Convergence Study 

In this section, first a convergence study for the above outlined 
method is presented. Then, the present method is illustrated for 
circular membranes with different loadings, shown in Fig. 3.6. 

1. Uniformly loaded circular membrane. 

2. A membrane with uniformly distributed load over the inner 
portion. 


(a) Uniform loading 



(b) Uniform load over inner portion 


n 



(c) Ring loading 


Fig. 3.6 Types of Loading on the Membrane 
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3. A membrane with a ring load. 

The stresses and displacements were expressed in dimensionless 
forms, using Hencky's normalized formulation as. 


771 


eq 


w 


w 


, ( Pa 1 1/3 
eq 


a 

r 


p 2 a 2 E 


171 


°e = 


(- 


°e 

22 ; i 7 i 
P a E 


eq 


) 


To study the convergence of the present method, the circular 
membrane was idealized into m number of regions with (m + 1) nodes, 
(where nodes are numbered from the center to the outside). The number 
of regions, m, used in this convergence study were 8, 16, 32, 60, and 
64. Figures 3.7(a) and 3.7(b) present the relative errors in the 
normalized maximum deflections and stresses, for a uniformly loaded 
membrane and for a membrane loaded uniformly over the region 
0 < — < 0.5. The solution shows rapid convergence and about 60 
regions were found to be necessary for a membrane loaded uniformly over 
the region 0 < ^ < 0.5, to yield a solution which is within 0.001 

percent of the converged solution. In contrast much fewer than 60 
regions were necessary for the uniformly loaded membrane. However, a 60 
region idealization is used in the analysis and all the results are 
presented for this idealization. 


The Newton-Raphson method used here needs initial 


a values for 
r 
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the iterative process. The sensitivity of the present method to these 

initial values was studied. Table 3.1 presents the number of iterations 

needed to achieve convergence for a range of initial values for the 

three configurations studied. The initial a values were as low as 

r 

three orders of magnitude smaller than the maximum converged value of 

a and as high as three orders of magnitude higher. For the uniformly 

loaded membrane with each of these initial values the present method 

converged to the same solution. The other two configurations showed 

similar convergence, thus showing insensitivity to the initial a, 

values. When the initial a values were farther away from the 

converged values the number of iterations needed were around 40 compared 

to about 8 when the initial cr^ values were closer to the converged 

a values, 
r 

3.4.1 Uniformly Loaded Circular Membrane 

The first problem analyzed was that of a uniformly loaded circular 
membrane, for which classical solution [29] exists. The uniform loading 
was represented by setting the magnitude of loading terms p^ through pgQ 
equal to unit values. Using these values in Eq. (3.22), the governing 
equation for the membrane was solved by using the procedure outlined 
earlier. The values of u and w displacements and stresses a and 
o Q are presented in Figs. 3.8 and 3.9, respectively. A comparison of 

U 

the present solution with Hencky's classical solution [29] and Kao and 
Perrone [53] nonlinear relaxation method solution is shown in Table 
3.2. The central displacement as well as central and edge radial 
stresses obtained by the present method are in excellent agreement with 
the earlier reported results as shown in Table 3.2. 
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3.4.2 A Membrane With Uniformly Distributed Load Over the Inner Portion 
The second example is that of a membrane with uniformly distributed 

load over an inner portion. A particular case of loading over the 

p 

region 0 < — < 0.1 was analyzed. In the corresponding solution, this 

loading was represented by setting the magnitude of loading terms p^ 

through pg equal to unit values and terms py through pg Q equal to 

zero. Again using the governing Eq. (3.22) with the loading terms 

discussed before, the solution was obtained for stresses and displace- 
ments at each node. The correspond! ng values of u and w displace- 
ments and o r and a Q stresses are shown in Figs. 3.10 and 3.11 
respecti vely . 

3.4.3 A Membrane With a Ring Load 

The last problem considered here is one where the membrane carried 
a uniformly distributed ring load. As a specific case the ring load was 
assumed to be spread over the region 0.5 < < 0.6. Since loading was 

a 

considered over the region 0.5 < < 0.6 in the corresponding 

a 

solution, this loading was represented by setting the magnitude of 
loading terms p 3 through p 3 g equal to zero, p 31 through p 3 g equal to 
unit values; and p 3 7 through pgg equal to zero. Using these values in 
Eq. (3.22), the governing equation for the membrane with the ring load 
was solved by using the procedure outlined earlier. 

The corresponding values of u and w displacements are presented 
in Fig. 3.12 and stresses o r and o Q are presented in Fig. 3.13. 
Figure 3.12 shows that the transverse displacement w is constant up to 
the ring load and then starts decreasing and becomes zero at the clamped 
edge. In contrast, the radial displacement u is zero at the center. 


Displacement w Displacement 




Fig. 3.10 Normalized Disj 
Loaded Over th 





Normalized stress, 
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increases gradually until it reaches the peak value in the loading 
region and then starts decreasing and becomes zero at the fixed edge. 
On the other hand, the normalized stresses and a Q are of the 

identical magnitude in the unloaded region 0 < — < 0.5 and then 

d 

decrease for larger values of r. 

3.5 Discussion 

As pointed out earlier, when inplane and transverse displacements 
differ by large amounts a simultaneous method of solution presents 
difficulties. As the present method does not use u and w as 

parameters, rather uses a single parameter these difficulties are 

avoided. Therefore, it is interesting to compare the differences in 

magnitudes of u and w values for various problems analyzed. For a 
membrane with uniformly distributed load u and w are of the same 

magnitude (Fig. 3.8), whereas in the case of membrane loaded over an 
inner portion and for the ring loaded membrane, the displacements differ 
by two (Fig. 3.10) and one (Fig. 3.12) order of magnitude, respectively. 
Because u and w do not differ by large amounts in the case of 

uniformly loaded membrane, one would expect the simultaneous solution to 
be efficient and feasible. Indeed, it is so as demonstrated by Kao and 
Perrone [53]. Because u and w differ by large amounts for the 

arbitrarily loaded membranes, one would expect difficulties with 
simultaneous solution method. This may be the reason for the limited 

numerical solutions for these two loading cases. The present method on 
the other hand avoided these problems by using a governing equation in a 
single parameter, the radial stress. The present method shows good 
convergence characteristics for all the problems studied and converges 
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to the correct solution even when the initial values are unrealistically 
large or small. Therefore, the nonlinear governing equation in terms of 
the radial stress and the Newton-Raphson technique appear to be the 
ideal choice for large deflection problems of arbitrarily loaded 
membranes. 


3.6 Concluding Remarks 

In this chapter the quasi-isotropic circular laminate with only 
extensional stiffness but no shear and flexural stiffnesses were modeled 
as circular isotropic membranes. A single nonlinear differential 
equation which governs the response of these circular clamped isotropic 
membranes under arbitrary axisymmetric loading was developed. This 
nonlinear equation was solved by using the finite difference method in 
conjunction with Newton-Raphson method. The numerical studies on the 
large deflection membrane analysis show that the present method of 
analysis yields accurate solutions for inplane and transverse 
deflections and stresses for the arbitrarily axi symmetri cally loaded 
circular membranes. 

In the next chapter a plate-membrane coupling model, which uses the 
small deflection plate solution obtained in Chap. 2 and the large 
deflection membrane solution obtained in this chapter, is presented to 
predict the large deformation behavior of thin circular quasi-isotropic 
laminates under point loads. 
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Table 3.1 - Sensitivity of the Present Method of the Initial 

Values of a “ 
r 


Initial 

a values b 
r 

Number of 

iterations required 

for convergence 

Uni form 

loading over 

the region 

0 < — < 
a 

1 

0 < ~ < 0.1 
a 

0.5 < I < 0.6 

a 

0.0001 

25 


22 

23 

0.001 

20 


16 

17 

0.01 

14 


10 

11 

0.1 

8 


8 

5 

1.0 

9 


19 

15 

10.0 

20 


27 

26 

100.0 

31 


41 

37 


Converged 0.4310 0.1104 0.1237 

a, (r =0) 


Converged 


0.0235 

DTiysn 

a (r = a) 
r 





a Membrane idealized with 60 regions 

^Constant radial stress a was assumed at all 61 nodes. 

r 
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Table 3.2 - Comparison of Normalized Displacements and Stresses 
for Uniformly Loaded Membrane 



Central transverse 
deflection 
w at r = 0 

Central radial 

stress 

a at r = 0 
r 

Edge Radial 
stress 

e r at r = a 

Hencky* 
(ref. 29) 

0.6536 

0.4310 

0.3280 

Kao and 
Perrone 
(ref. 53) 

0.6541 

0.4289 

0.3306 

Present 
Resul ts 

0.6534 

0.4310 

0.3329 

*Values taken 

from reference 53. 





Chapter 4 


PLATE-MEMBRANE COUPLING MODEL 
4.1 Introduction 

In Chapters 2 and 3, numerical solution methods to analyze the 
small deformation behavior of circular quasi-isotropic plates and the 
large deflection behavior of circular isotropic membranes are 
described. In this chapter, a plate-membrane coupling model to study 
the large deformation behavior of clamped circular plates is 
developed. The model developed herein uses the plate and membrane 
solutions in conjunction with a coupling principle. A numerical 
solution method to analyze the clamped circular quasi-isotropic 
laminates under point load, using pi ate -membrane coupling model is 
presented. 


4.2 Analysis 

In the plate-membrane coupling model two different plate problems 
are analyzed. First, a thin plate with shear and flexural stiffnesses 
but no mid-plane extensional stiffness under point load is considered. 
The deflected shape of the plate is then obtained by using the classical 
small deformation theory. Second, a plate with mid-plane extensional 
stiffness, but no shear and flexural stiffnesses, i.e., a membrane is 
considered. The deformed plate shape obtained in the plate problem is 
used in the membrane problem. For this deformed shape membrane loads 
are calculated by using a nonlinear membrane theory. Large deformation 
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solutions of thin plates are then obtained by coupling these solution of 
plate and membrane problems. The details of the plate-membrane coupling 
model are given below: 

4.2.1 Plate Configuration 

Consider a clamped circular plate of thickness h and radius a, 
subjected to a central point load P as shown in Fig. 4.1. The plate is 
a quasi-isotropic laminate of T300/5208 graphite/epoxy material with 
stacking sequence [45/0/-45/90] s and material properties given in Table 
2 . 1 . 

To analyze this quasi-isotropic clamped circular plate for large 
, deformation behavior, first a circular plate with a prescribed central 
deflection w 0 is considered. The objective of the analysis is to 
determine the large deformation shape of the plate w(r) and the 
central concentrated load P. 

This problem as mentioned earlier, can be decomposed into two 
component problems: problem-1 is a plate with shear and flexural 

stiffnesses but no mid-plane extensional stiffness, problem-2 is a plate 
with mid-plane extensional stiffness but no shear and flexural 
stiffnesses, i.e., a membrane. 

4.2.2 Problem-1: Plate Problem 

In problem-1 the plate was assumed to be flexurally isotropic even 
though a quasi-isotropic plate is anisotropic in flexure. (This 
assumption was experimentally verified and will be discussed in detail 
in Chap. 6.) 


From Chap. 2, the governing differential equation for axisym- 
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metrical bending of a clamped circular plate subjected to arbitrary 
axisymmetric loading and undergoing small deformation is (Eq. (2.12)-) 


d 3 w + JL d 2 w 1 dw Q(r ) 
dr 3 r ^ ° D e, 


(4.1) 


where Q(r) is a shear force at any radius r and is given by (Eq. (2.13)) 


2 n r Q(r ) = 2 w l p(0 d£ 


(4.2) 


where p(£) is the intensity of loading at any radius £. 

In equation (4.1), D g q is an equivalent flexural modulus for the 
quasi-isotropic laminate. As discussed in Chap. 2, the equivalent 

modulus D e q was obtained by equating bending energies of the quasi- 
isotropic laminates and an equivalent isotropic plate. Equation (4.1) 
can be solved to obtain the deformed shape of the plate for a given 
intensity of load p(r). 


4.2.3 Problem-2: Membrane Problem 

In problem-2, a plate with mid-plane extensional stiffness but no 
shear and flexural stiffness, i.e., a membrane is analyzed. Since the 
quasi-isotropic laminates are axially isotropic, the governing 
differential equations for these laminates are essentially the same as 
that for a circular isotropic membrane. From Chap. 3, the governing 
equations for the large deflections of circular isotropic membranes 
subjected to arbitrary axisymmetric loading are (Eqs. (3.18) and (3.19)) 



do 

cfr~ 


+ r 


d 2 o. 


dr‘ 


r ] + (£) 2 


0 


(4.3) 


and 
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f 

fdw>|2 r'o 
l 3r J " L ~ 


p(S) S d£ 2 

r h a ^ 
r 


(4.4) 


where 0 p is the radial stress. If the displacement w is known, the 
rotation (^-) is known and hence the radial stress can be determined 
by integrating Eq. (4.3). With this radial stress a distribution, 
the loading p{r) on the membrane can be determined by using equation 
(4.4). 

4.2.4 Coupling of the Two Problems 

The large deformation solution for quasi-isotropic clamped circular 
laminate with the prescribed central point deflection w Q can be obtained 
by using the governing equations, Eqs. ( ( 4 . 1 ) -( 4 . 4 ) ) and by the 
appropriate coupling of plate and membrane problems as shown in a flow 
chart of the Fig. 4.2. This procedure is outlined in Fig. 4.3, and is 
as follows: 

The central point load P and the deformed shape w(r) for a 
clamped circular plate were obtained by using the prescribed central 
deflection w Q and the governing equation, Eq. (4.1) (steps 1 and 2 in 
Figs. 4.2 and 4.3). This deformed shape w(r) was used in the membrane 
problem, to determine the membrane loads p(r) (step 3 in Figs. 4.2 and 
4.3). The loads which are equal in magnitude to p(r) but opposite in 
direction, and the point load P were applied to the plate problem (step 
4 in Figs. 4.2 and 4.3). For this new loading the new deformed shape 
w' (r) of the clamped circular plate was obtained by solving the plate 
problem (step 5 in Figs. 4.2 and 4.3). Since the membrane loads p(r) 
were applied in the opposite direction, this caused reduction of the 
central deflection of the plate by an amount A w Q (step 6 in Fig. 4.2). 
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Iterate 
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Fig. 4 


STEP 1: 


I-*- STEP 2: 


STEP 3: 


STEP 4: 


STEP 5: 


1 STEP 6: 


STEP 7: 


STEP 8: 


*• STEP 9: 



3 Flow Chart for the Plate-Membrane Coupling Model 







86 


An incremental load AP is necessary to nullify the reduction of 
the central deflection, A w Q . The magnitude of the incremental load 
AP, is calculated by using Eq. (4.1) (step 8 of Fig. 4.3). This 
incremental load AP, together with the original point load P and the 
membrane loads which are equal in magnitude to p(r) but opposite in the 
direction were applied to the plate problem (step 9 of Figs. 4.2 and 
4.3). This load system yields a central deflection equal to w Q , the 
original value. However, the new deflected shape w(r), determined from 
Eq. (4.1) is in general different from the original deflection shape. 
This iterative procedure is repeated until the deflected shapes obtained 
in any two consecutive iterations are almost identical, i.e. step 7 in 
Fig. 4.3 is satisfied. 

A complete solution can now be obtained by coupling the plate- 
membrane solutions. The coupled plate and the membrane are analogous to 
two coupled parallel springs, one with linear (plate) stiffness and 
other with a nonlinear (membrane) stiffness. If these springs undergo 
the same deflection w Q then the unknown load P applied to this spring 
system is the sum of the loads carried in the two springs. Similarly in 
the present coupling model the deflected shapes of the plate and the 
membrane are identical and therefore the loads are additive. 

V/hen the deflected shapes obtained in any two consecutive 
iterations are nearly identical the plate carries a central point load 
and a distributive interactive load ( — p { r ) ) and the matched membrane 
carries the opposite of distributive interactive loads p(r). When the 
plate and membrane are coupled, interactive loads on the plate and 
membrane nullify each other and the total load is simply a single 
central point load on the plate. The complete plate-membrane solution 
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thus yields the large deflection shape of the plate and the magnitude of 
the central . concentrated load, for a prescribed central deflection w Q . 

Instead of prescribing a central deflection w Q , one could formulate 
the problem for a prescribed central point load P. The analysis 
procedure for this case is straight forward with minor differences and 
is as follows: 

The problem as before, is also decomposed into two component 
problems. By using the center point load P in the plate problem, 
deflected shape w(r) and the central deflection w Q is obtained by using 
Eq. (4.1). Then steps 2 through 6 (Figs. 4.2 and 4.3) are repeated as 
done earlier. Since in step 4, the membrane load p(r) is applied in the 
opposite direction, this causes a reduction in central deflection by 
A w Q (step 6 of Fig. 4.3). At this stage instead of incrementing the 
center point load P, the new membrane loads p(r) are calculated by 
using the deformed shape w'(r) obtained in step 5 (step 3 of Fig. 4.3). 
This procedure is repeated until the deflected shapes obtained in any 
two consecutive iterations are almost identical. 

A complete solution can now be obtained by coupling the plate- 
membrane solutions. When the final plate and membrane solutions are 
coupled as before, the interactive loads on the plate and membrane 
nullify each other and the original prescribed central point load is 
left on the plate. The corresponding matched central deflection w Q , and 
the deflection shape w(r), are the large deformation solutions of the 
clamped circular plate subjected to a prescribed central point load. 

The above procedure, although presented for a central concentrated 
load, can be applied in a similar manner to obtain the large deformation 
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solutions for arbitrarily axisymmetrically loaded clamped circular 
plates. If, on the other hand, the large deformation shape of the plate 
is known, it can be used in the linear plate theory to calculate the 
plate loads, and in the nonlinear membrane theory to calculate the 
membrane loads. The sum of these two loads gives the complete large 
deformation solution. 


4.3 Solution Method 

The plate governing equation, Eq. (4.1) for the deflecton w, was 
solved by using finite difference method as described in Chap. 2. The 
ordinary differential Eq. (4.1), was replaced by a set of linear 
algebraic equations using finite difference quantities. These linear 
algebraic equations were solved numerically, to obtain rotations and 
displacements. The nonlinear governing equations, Eqs. (4.2) and (4.3) 
were replaced by a set of nonlinear algebraic equations using finite 
difference quantities. These nonlinear equations were solved by using a 
Newton-Raphson method in conjunction with the finite difference method, 
as described in Chap. 3. 

Using these plate and membrane solutions in conjunction with the 
plate-membrane coupling model discussed before, the large deformation 
shapes for the clamped circular quasi-isotropic laminates were obtained. 

From the plate-membrane coupling model it was observed that as 

plate central deflection to plate thickness ratio increases, the number 

of iterations required to obtain the large deformation solution 

increases. Table 4.1 presents the number of iterations required to 

obtain the large deformation solution of the clamped circular quasi- 

w 

isotropic laminates under point loads for (— pj-) ratios ranging from 0.5 
to 2.0. 
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4.4 Results and Discussion 

In this section analytical deflected shapes obtained with the 

plate-membrane coupling model, for the central deflection to plate 

w 

thickness ratios ranging from 0.5 to 2.0 are presented. These 

analytical deflected shapes obtained by using the plate- membrane 

coupling model are compared with the deflected shapes obtained by using 

the classical solution [27]. Load-displacement curve obtained by using 

the plate-membrane coupling model is compared with the classical 

solution [27]. Figures 4.4 and 4.5 show the analytical deflected shapes 

obtained by using the plate-membrane coupling model; for (— pj-) ratios 

of 0.5 to 1.0, and 1.5 and 2.0. Figures 4.4 and 4.5 indicate that the 

functional form of the deformed shapes of the plate is different for 
w o 

various (— ^) ratios. 

To compare the plate-membrane coupling analysis, the classical 
solution [27] based on the energy method was considered. Appendix B 
presents the large deformation classical solution based on the energy 
method for a clamped circular plate subjected to the center point load 
P. The classical solution assumes that the functional form of the 
deformed shape of the plate is identical to the functional form of the 
deformed shape of the plate determined by the small deflection 
solution. Thus as per the classical solution [27], the large deformed 
shape of a clamped circular plate under center point load is given as: 

w(D = « [1 - (l!) + 2 An (I)] (4.5) 

a a 

The transverse displacements w obtained by using the plate- 
membrane coupling model and the classical solution (Eq. (4.5)) for 
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various (— ^) ratios ranging from 0.5 to 2.0 are presented in Table 

4.2. The results obtained by using the plate-membrane coupling model 

w 

show that the functional form of the deflected shapes for various (— ^-) 
ratios is different. In contrast, as mentioned earlier the classical 
solution [27] based on energy method, assumes that the functional form 
of the deflected shape of the plate is identical for all (— ratios. 
The deflected shapes obtained by using the plate-membrane coupling model 
are compared with those obtained by using the classical solution [27]. 

Figure 4.6 presents the percentage errors, in the deflected shape 

of the classical solution relative to the present solution, at various 
w o 

radii for (— ^-) ratio of 2.0. It is observed that the maximum error is 

about 4 percent and occurs at (-^-) = 0.233. Although relative errors 

in the displacements are small, these errors get magnified if the 

curvatures obtained by the plate-membrane coupling model and the 

classical solution are compared. Figure 4.7 compares the curvatures 

obtained with the plate-membrane coupling model and the classical 
w 

solution, for (— pp) ratio of 2.0. Figure 4.7 shows that the plate- 
membrane coupling model predicts higher values of curvatures than those 
obtained by using the classical solution up to (— ) ratio of 0.4 and 

a 

beyond (— ) = 0.9. However, in the range 0.4 < ^ < 0.9, the 

classical solution predicts higher values of curvatures. Since stresses 
are proportional to the curvatures even small differences in the 
curvature result in significant errors in stress predictions in the 
plate undergoing large deformation. 

In the plate-membrane coupling model, another important point that 
was noticed was the movement of the radius of inflection (radius at 

which curvature is zero), as a function of the central deflection to the 
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w 

plate thickness ratios (-^-) . The classical solution [27], assumes 

that the functional form of the deformed shape is identical for all 
w 

ratios, hence the radius of inflection shows no movement with the 
change in ratios and is located at (•£■) = 0.3678. (The detail 

derivation for the radius of inflection is presented in Appendix C.) 
The plate-membrane coupling model on the other hand does not make any 
assumption on the deformed shape of the plate and hence on the radius of 
inflection. The deformed shape and the radius of inflection are 
obtained as part of the solutions. Figure 4.8 shows the movement of the 
radius of inflection with the change in the central deflection to the 

W Q 

plate thickness ratios (— jr) . The symbols in Fig. 4.8 show the 

w 

position of the point of inflection for different values of (— ^-) . For 

increasing central deflection the point of inflection moves toward the 

center of the plate (r =0). The inward movement of the point of 

inflection with increasing (— ^-) ratio occurs, since, as ratio 

increases, the plate behaves more like a membrane than a flexural plate, 
w 

That is as (— ratio increases, the external load is equi 1 ibriated 
more and more in membrane action than in flexural action. 

Figure 4.9 presents a comparison between the load-displacement 

curves obtained by using the plate-membrane coupling model and the 

classical solution [27]. This figure indicates that the classical 

solution based on the energy method predicts the load-displacement 

w 

behavior accurately up to (— ^) ratio of 1.0. However for higher 
w 

values of (— ^) ratios the classical solution yields lower 

displacements than the plate-membrane coupling solution. 



.4 


Classical solution 


Fig . 4 



w 0 /h 


8 Movement of Radius of Inflection as a Function of 
Center Deflection 



0 5 10 15 20 25 30 


Fig. 4.9 Load-Deflection Curve Comparison 
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4.5 Concluding Remarks 

In this chapter circular quasi-isotropic composite plates were 
analyzed to study the large deformation behavior under point loads. A 
pi ate -membrane coupling model was formulated and solved numerically to 
obtain the large deformation shapes for thin circular composite 
laminates under the point loads. These large~deformation shapes were 
compared with the large deformation shapes obtained by using the 
classical solution [27] which assumes that the functional form of the 
deformed shape of the plate is identical to the functional form of the 
deformed shape of the plate determined by the small deflection plate 
solution. The comparison shows that the functional form of the deformed 
shape of the plate undergoing large deformations is different from the 
small deflecton plate solution in that the deformed shape is a function 
of the center point displacements and thus is different for different 
load levels. The classical solution and the plate membrane coupling 
solution are in good agreement up to the central deflection-to-plate 
thickness ratio of 1.0. For higher values of the center deflection-to- 
plate thickness ratios the classical solution yields lower displacements 
than the plate-membrane coupling solution. This plate membrane coupling 
model is experimentally verified in Chap. 6. 
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Table 4.1 - Number of Iterations to Obtain Large Deflection 

Solution of a Clamped Circular Plate Under Central 
Concentrated Load Using the Plate-Membrane Coupling 
Model . 


Plate center deflection 

w 

to thickness ratio (— ^-) 


Number of 
i terations 


0.5 

9 

1.0 

11 

1.5 

14 


2.0 


21 




TOO 


Table 4.2(a) - Comparison of the Transverse Displacements w 

Obtained by Using the Pl^te-Membrane Model and the 

Classical Solution for Ratios of O.S and 1.0. 


r/a 


Classical 

Solution 


Plate-Membrane 
Coupling Model 




0.0 

0.50000 

0.50000 

0.1 

0.47197 

0.46764 

0.2 

0.41562 

0.40842 

0.3 

0.34664 

0.33791 

0.4 

0.27339 

0.26484 

0.5 

0.20171 

0.20311 

0.6 

0.13610 

0.13506 

0.7 

0.08023 

0.07762 

0.8 

0.03719 

0.03828 

0.9 

0.00966 

0.01001 

1.0 

0.00000 

0.00000 

0.0 

1.00000 

1.00000 

0.1 

0.94394 

0.93050 

0.2 

0.83124 

0.80831 

0.3 

0.69328 

0.66745 

0.4 

0.54678 

0.52491 

0.5 

0.40342 

0.39392 

0.6 

0.27220 

0.27028 

0.7 

0.16046 

0.16570 

0.8 

0.07438 

0.07719 

0.9 

0.01932 

0.02074 

1.0 

0.00000 

0.00000 
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Table 4.2(b) - Comparison of the Transverse Displacements w 

Obtained by Using the Pl^te-Membrane Model and the 

Classical Solution for Ratios of 1.5 and 2.0. 


r/a Classical 

Solution 


PI ate -Membrane 
Coupling Model 




0.0 

1.50000 

1.50000 

0.1 

1.41592 

1.38595 

0.2 

1.24686 

1.19548 

0.3 

1.03992 

0.98512 

0.4 

0.82018 

0.78421 

0.5 

0.60513 

0.58841 

0.6 

0.40830 

0.40245 

0.7 

0.24068 

0.25309 

0.8 

0.11156 

0.12140 

0.9 

0.02897 

0.03286 

1.0 

0.00000 

0.00000 

0.0 

2.00000 

2.00000 

0.1 

1.88789 

1.83399 

0.2 

1.66248 

1.57111 

0.3 

1.38656 

1.29297 

0.4 

1.09357 

1.02026 

0.5 

0.80685 

0.77129 

0.6 

0.54441 

0.53813 

0.7 

0.32091 

0.34512 

0.8 

0.14875 

0.16870 

0.9 

0.03863 

0.C4561 

1.0 

0.00000 

0.00000 



Chapter 5 


FINITE ELEMENT MODEL 
5.1 Introduction 

In the previous chapter, a plate-membrane coupling model was 
presented to obtain the large deformation behavior of circular quasi- 
isotropic laminates under point loads. To verify the plate-membrane 
coupling model, a finite element formulation, which uses appropriate 
linearization of strain-displacement relations [36] is developed and 
presented in this chapter. The formulation is based on a variational 
technique and uses the principle of minimum potential energy [28]. As 
pointed out in Chap. 1, several investigators [36-44] used a lineariza- 
tion technique, to study the geometrically nonlinear circular plate 
problems. However, earlier investigators ignored the membrane effects 
due to mid plane stretching. In this chapter, the membrane effects due 
to the mid-plane stretching are incorporated in the finite element 
model. This model is used to analyze the large deformation behavior of 
the circular quasi-isotropic laminates under point loads. The results 
from this model are compared with the results from the plate-membrane 
coupling model. 


5.2 Finite Element Method 

The basis of the finite element is the representation of a body or 
structure by an assemblage of subdivisions (finite elements). Simple 
functions are then chosen to approximate the distribution or variation 
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of the actual displacements over each element. These functions are 
usually referred .to as shape functions. A variational principle, such 
as the principle of minimum potential energy, is then employed to obtain 
the set of equilibrium equations for each element. The equilibrium 
equations for the entire body are then obtained by combining the 
equations of the individual elements. The equations are modified for 
the given force or displacement boundary conditions and then solved to 
obtain the unknown displacements. 

Mathematically, the finite element representation of continuum 
strongly resembles the Ritz method [ 28 ], in which the displacements of 
the plate are approximated by the sum of the functions, each multiplied 
by an unknown constant. These unknowns are determined from the minimum 
potential energy theorem. While using the Ritz method, the assumed 
series expression describes the total displacement field of the entire 
plate, in the finite element method individual displacement patterns for 
each element are assumed. The entire displacement field of the plate 
can be approximated piecewise. The total potential of the plate, 
obtained from the sum of the total potentials of the individual 
elements, has a stationary value when the node points are in 
equilibrium. This conditions leads to minimization of the total 
potential of the structural system, which, in turn, yields the 
displacement field corresponding to its equilibrium condition. 

The application of the finite element method to obtain the large 
deformation solution for the quasi-isotropic circular plates under point 
loads is presented in following sections. 
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5.3 Plate Configuration 

As before the plate configuration used was that of a quasi- 
isotropic circular plate with a concentrated load at the center. The 
plate is assumed flexurally isotropic and the equivalent flexural 
modulus D e q of the plate is computed as described in Chap. 2. With this 
assumption of the flexural isotropy the problem of a circular quasi- 
isotropic plate subjected to a central point load reduces to an 
axisymmetric problem. This problem is then analyzed by the finite 
element method as shown in the following sections. 


5.3.1 Strain-Displacement Relations 

For the axi symmetrically loaded circular plate the deflection 
surface is axisymmetrical . Due to the symmetry the shearing stresses 
x r0 are zero. Also from the basic assumptions of the plate theory 
(Chap. 2), shearing stresses x and normal stresses o ■ are zero. 

For the circular plate undergoing large deformations, the radial 
and tangential bending strains at any distance z from the middle surface 
are, 

= z 
bendi ng 




d w 

\ 




F? 

1 


e 

r 1 

1 

dw 

{ = 

z 

e 0 j 

r 

W 

I 



and the membrane strain-displacement matrix is (Eqs. (3.3) and (3.4)) 


membrane 


du 1 rdw}2 

Fr 7 l 3r J 

u 
F 


| 

e 


r 

1 


1 

e 


(5.2) 


m 


Therefore the total radial and tangential strain-displacement 
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matrix is 


e 


e 



e 

r 

= Z 

r 


+ 

r 

CD 

W 


e 6 

b 


E e 


or in concise form 

{e} = z (e. } + (e } 
b m 


( 5 . 3 ) 


( 5 . 4 ) 


The stresses and strains are related by a material properties 
matrix, [H], forming the following constitutive equation 


{a} = [H] (e) 


(5.5) 


where 


[H] 


"eq 


(1 ~W 



(5.6) 


where E g q and v are the equivalent Young's modulus and Poisson’s 
ratio for the quasi -i sotropi c plate and were calculated by using the 
technique presented in Chap. 2. 


5.4 The Axi symmetric Finite Element 

An axisymmetric finite element in the form of a ring of constant 
cross section is shown in Fig. 5.1. The node points of such an element 
are in fact nodal circles, and the volume of such an element is 
dependent on both its cross-sectional area and the radii of these nodal 


circles. 
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A convenient approach for derivation of the finite element 

governing expressions is based on the principle of minimum potential 
energy. Consider a circular plate discretized into n number of 

axi symmetric elements. Each element is of a constant thickness h, with 
r^ and r 0 as inner and outer radii (Fig. 5.1). Let U e be the strain 

energy of the axisymmetric element and A represents the surface area 

of the element. The total potential energy Tc g of an element is the 
total strain energy of the element minus the total work done and can be 
expressed as 

* - U - // (S) T {f> dA (5.7) 

e e A 

where { S) is a vector of nodal displacement and {f } is the 

corresponding load vector containing inplane loads I, the transverse 

loads p and the applied moments M at all nodes on the element. U e in 

the Eq. (5.7) is the strain energy of the axisymmetric element and is 

given by 

U = -i / {a}^ {e} d (vol ) (5.8) 

e 6 vol 

vol in Eq. (5.8) is the volume of the axisymmetric element. 
Substituting Eq. (5.5) into Eq. (5.8) strain energy for the element was 
written as 

U i J (e} T [H] (e) d (vol) 
e £ vol 


(5.9) 
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where [ H 3 is the constitutive matrix given by Eq. (5.6). For a plate 
under combined bending and stretching the strain-displacement relation 
is given by Eq. (5.4). Substituting Eq. (5.4) into equation (5.9), 
strain energy U e is: 


U 


e 


n C 


eq 



V T [Q] {c m } r dr 


(5.10) 


+ it 




{e b > T [Q] {e b > r dr 


C eq and D eCj in Eq. (5.10) are called the membrane and the bending 
stiffness coefficients respectively and are given by 


and 


eq 


E e, h 

(1-V 2 ) 

eq 


and 


£ h' 
eq 


eq 


12 (Kq* 


(5.11) 


CQ3 


1 

V 

eq 


v 

eq 

1 


(5.12) 


Substituting Eq. (5.10) into Eq. (5.7), the total potential energy 
of the element is: 


"e = C eq <E m ,T tQ] U m > r dr 
+ it / r ° {t^} 1 [03 U b ) r dr 

- 2u J r ° (S) T {f ) r dr] 


(5.13) 
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By using the principle of minimum potential energy, the element 
stiffness matrix was obtained as follows: 

5.5. Stiffness Matrix of an Axisymmetric Element 

A cubic displacement in w and a linear displacement in u were 
assumed in the plate element as 

2 3 

w = a 0 + a^r + a 2 r + a 3 r (5.14) 

u = b Q + b.r (5.15) 

These displacement functions ensure the continuity of displacements 
between adjacent elements. Thus these displacement functions satisfy 
the compatability requirements, which state that at element interfaces 
the field variables u, w and any of its partial derivatives up to one 
order less than the highest order derivative appearing in the energy 
expressions itp must be continuous. 

In the plate problems, the bending strains are defined by second 
derivatives of the transverse displacements. Therefore a cubic function 
in w will give a constant strain in an element. Similarly the membrane 
strains are defined by first derivatives of the radial displacements and 
thus a linear function in u will give a constant strain in an element. 
Therefore a cubic displacement funciton in w and a linear displacement 
function in u meets the constant strain condition in the element and 
thus satisfies the convergence criteria in the finite element. 
Equations (5.14) and (5.15) contains six unknowns a Q , a^, a 2 , * 3 , b Q , 
and b^ and the axisymmetric plate element has two nodes, therefore each 



no 


node was assumed to have three degrees of freedom, the radial displace- 

dw 

ment u, the transverse displacement w, and the rotation . These 3 
degrees of freedom satisfy the compatability requirements discussed 
earlier, which requires that the u, w, and should be continuous at 

element interfaces. 

For any axi symmetric element with radius r^ and r 0 , nodal displace- 
ments can be expressed in terms of generalized displacements as: 



(5.16) 


where subscripts i and o represent inner and outer nodes. Equation 
(5.16) was expressed in the concise form as follows: 


{ 6> e = [A] {a} 


(5.17) 


where (6> e = the nodal displacement vector. From the foregoing, the 
solution for the unknown constant is 


(a) = [A]" 1 {6} 

e 


(5.18) 


Now consider the total potential energy equation (Eq. (5.13)). Applying 
the principle of the minimum potential energy, the variation in the 
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total potential energy should be equal to zero. Therefore the variation 
in the potential energy from Eq. (5.13) is 


A* e - [2 * D, q / r ° U Cj,) 7 [Q] {£ b ) r dr 


+ 2n C / 0 (Ae } T [Q] {e > r dr 
eq r.. m m 


2n / 0 {As {f } r dr] = 0 
r i 


(5.19) 


The bending strains (e^) in Eq. (5.19) are (Eq. (5.1)) 


V = 


2 

d^w 

— ? 
dr 

_1 dw 
r dr 


Using the shape functions (Eq. (5.14)), {^] and can be 

u ‘ dr 


2 

-d^w- 


written as 


and 


dw 

W ' a l 


+ 2 a 2 r + 3 a 3 r 


2 



2 a£ + 6 a^r 


(5.20) 


(5.21) 


Using Eqs. (5.20) and (5.21), (e^) can be expressed in matrix form as 
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{e b } = 


0 

0 


or concisely as, 


0 

-1 

r 


-2 -6 r 0 0 

-2 -3 r 0 0 





J 


(5.22) 


U b ) = [R] {a} 

= [R] [A]' 1 {6} 

e 

Using Eq. (5.24), (Ae^ can be written as 
{Ae b > = [R] [A]" 1 {A6) e 


(5.23) 

(5.24) 

(5.25) 


5.5.1 Nonlinear Terms: - Linearization Procedure 

Consider the membrane strains {e^} in the equation (5.19). They 

are: 

dr 1 rd w->2 

Hr 7 l dr J 

£ 

2 

Since the vector {e } contains a nonlinear term -i- (4^-)^ a 

m l v dr 

linearization technique [36] needs to be used. The linearization 

procedure used is as follows: 

Consider a nonlinear term in the membrane strains {e m > . 



Let 


1 dw 
7 Hr ~ 


(5.26) 
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where B is called as a linearization constant. Using this definition 
the membrane strains can be written as. 


UJ = 
m 


dr o dw 
W + 
u 
r 


Using Eq. (5.15), and can be written as 


(5.27) 


du K 
W b l 


(5.28) 


and 


“=_° + b. 
r r 1 


(5.29) 


Using Eqs. (5.20), (5.28) and (5.29), (e m } can be expressed in matrix 
form as 

r 

a o 
a l 


UJ = 

m 


0 B 2 Br 3 Br 
0 0 0 0 


0 

1 




(5.30) 


or concisely as. 


{e } = CT] 

m 


(5.31) 


= [T3 [A]" 1 (6) 


(5.32) 
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The next step is to evaluate the variation of the membrane strains, 

{Ae } . This is achieved as follows: 
m 

The membrane strains {e } can be represented as the sum of the two 

m r 

vectors, one containing only the linear terms, whereas the other 
containing the nonlinear terms. Therefore, 


where 


(e } = {e, } + (e... } 
m L NL 


{e^} 


du 

Hr 

F 


{ v 


1 r dw 'i 2 

7 > 


The variation of {e }, then is 

m 


{Ae } = {Ae } + {Ae... } 
m L NL 


( 5 . 33 ) 


( 5 . 34 ) 


( 5 . 35 ) 


( 5 . 36 ) 


or 


or 



( 5 . 38 ) 
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The linearization constant B is then used in Eq. (5.38) and (Ae } is 

^ m 

rewritten as. 


(Ac } = 

m 



(5.39) 


Using the shape funcitons for u and w the variation of membrane strains 
can be expressed in the matrix form as. 


r 



+ 



2 B 4 Br 
0 0 


2 

6 Br 
0 


0 

0 


0- 

0 - 


r 



} 

(5.40) 


or in concise form as. 


{Ae m > = A [[Y 1 ] [A]" 1 {6) e + [Y^ [A]" 1 (6> e J 


(5.41) 


= [Y 1 + Y 2 ] [A]" 1 (6) e 


(5.42) 


Lastly consider a generalized displacement vector { S > 
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(5.43) 


Using Eqs. (5.14), (5.15) and (5.20), a generalized displacement vector 
can be expressed in the matrix form as 


(S) = 


0 0 

1 r 

0 1 


0 0 1 r 

2 3 

r^ r J 0 0 

2r 3r 2 0 0 



or in concise form 


(5.44) 


(S) = [N] {a} 


(5.45) 


By using Eq. (5.18) for {a}, in Eq. (5.45), Eq. (5.45) was written as 

(S> = [N] [A]" 1 {6} (5.46) 

e 

and 

{AS} = [N] [A]' 1 (A6> e (5.47) 

In summary the matrices obtained so far are. 


U. } = [R] [A]' 1 { 6) (Ac.) = [R] [A]" 1 (A6) 

d e d e 


UJ = [T] [A]" 1 { 6> e {Ae m } = [Y 1 + Y ? ] [A]" 1 {A6> e 
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CS> = [N] [A]' 1 { 6 } e {AS} = [N] [A]' 1 {A6) e 

The variation in the potential energy of the element (A ) of 
Eq . (5.19) is then, 

r 

A * e = 271 D eq ^r° {A6} J [R]T ^ {6 > e r dr 

p 

+ 271 C eq ^r° {A6} I ^ a3_1 ] T ^ Y i + Y 2 ]T CQ] [T] [A] ' X {6} e r dr 

p 

- 2n / r ° {A6}^ [CA]" 1 ] 1 [N] T {f > r dr = 0 (5.48) 

Since Eq. (5.48) should be valid for arbitrary values of {A6) e> Eq. 
(5.48) reduces to 

([k L ] e + ^NlV (6) e = {F, e (5 - 49) 

where [k^] e and ^NL^e are ca ^ ec ^ the bending and the membrane stiffness 
matrix for the element and {F) g is the corresponding load vector and 
are given by 

p 

[k L ] e = 2u D eq ^r° U A ]j ^ ^ t A ] _1 r dr (5.50) 

Ck NL ] e = 2 * C eq 4° tCA3_1 [Y 1 + Y 2 ]T CQ] [T] CA] " 1 r dr 


(5.51) 
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p 

(F> = 2n J ° [CAD" 1 ] T [N] T {f} r dr (5.52) 

e r. 

Assembling the element stiffnesses and nodal loads, global stiffness 
matrix and total load vector can be obtained and Eq. (5.49) can be 
written as: 


[K] {6} = (F) 


(5.53) 


where [K] is a global stiffness matrix and is given by 


[K] = I 
1 


«iA + 


t k NL V 


(5.54) 


and 


n 

{F} = E {F} (5.55) 

1 e 

5.5.2 Boundary Conditions 

For the clamped circular quasi-isotropic plate with axisymmetric 
loading, boundary conditions are 

( dVJ \ 

ppj equal zero at 

the center (r = 0). 

2. The radial displacement u, the transverse displacement w, and 

( dv/ ^ 

-gpj equal zero at the clamped edge (r = a). 
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5.6 Solution Procedure 

To solve the governing Eq.' (5.53) for the nodal displacements 

{6}, the solution domain was discretized into n axisymmetric elements 

and (n+1) nodes. Each node had three degrees of freedom, the radial 

f dw ^ 

displacement u, the transverse displacement w, and the rotation . 

In Eq. (5.53), [K] is a global stiffness matrix and is given by the Eq. 
(5.54) as: 


CK] = Z {[k. L + 
1 L e 


K«.V 


In order to evaluate the global stiffness matrix [K], it is 
required to obtain the bending stiffness matrix [k|_] e and the membrane 
stiffness [k NL ] e f° r each element. However, the membrane stiffness 
matrix [k NL ] e ^9* (5.51)) contains the linearizing function "B" and as 
the function "B" was not known a priori, an iterative scheme was adopted 
in the present solution. The details of the iterative scheme are given 
below and presented in the flow chart of Fig. 5.2. 

The linear bending stiffness [k|_] e and the load vector (F) e were 
computed by using Eqs. (5.50) and (5.52) (step 3, Fig. 5.2). Since the 
linearizing function B, was not known a priori, it was assumed zero and 
hence the membrane stiffness [ kfj|_ was a null matrix (step 4, Fig. 
5.2). Using Eqs. (5.54) and (5.55), the global stiffness matrix [K] and 
the load vector { F} were computed (step 6, Fig. 5.2). Equation 
(5.53) was then solved to obtain the linear displacement solution {6} 
(step 7, Fig. 5.2). By using the linear displacements and the rotations 
at each node, the linearizing function "B" and the membrane stiffness 
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Fig. 5.2 


Flow Chart for the Finite Element Model 
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[k N j_] e were computed (step 10, Fig. 5.2). The total stiffness matrix of 
the element [k] e is then the sum of the linear bending stiffness [kj_] e 
and current membrane stiffness [k^^. The global stiffness matrix [K] 
was then formed and the new displacements were computed. Using the new 
displacements a new [k N |_] e matrix is computed. This iterative procedure 
was continued until there was no significant difference in displacements 
between the successive iterations. 

5.7 Convergence Study 

To study the convergence of the present method, the circular plate 
was discretized into n axisymmetric elements and (n+1) nodes. A 
systematic convergence study was made by increasing the number of 
elements in the idealizations. The number of elements used in this 
convergence study were 2, 4, 8, 16, 32, 40, and 48. For the centrally 
loaded clamped circular plate with center deflection to plate thickness 
ratio of 2.0, a 48 element idealization was found to obtain a converged 
nonlinear solution. Figure 5.3 presents the relative errors in the 
maximum deflections normalized with respect to the converged maximum 
center deflections. The solution shows a rapid convergence and about 32 
regions were found to be necessary to yield a solution which is with 
0.01 percent of the converged solution. 

By using the 48 element idealization, the clamped quasi-isotropic 
circular plate was analyzed to obtain the large deformation solutions 
for various center def lection-to-plate thickness ratios. As the center 
deflection-to-plate thickness ratio increases, the number of iterations 
required to obtain the large deformation solution increased. Table 5.1 
presents the number of iterations required to obtain the large deforma- 
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Fig. 5.3 Convergence Study for the Plate Carrying a 
Concentrated Load at the Center 
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tion solution of the clamped circular quasi-i sotropic laminates under 
point loads for central defied ton-to-plate . thickness ratios ranging 
from 0.5 to 2.0. 


5.8 Results and Discussion 

In this section, the deflected shapes obtained with the finite 
element model for the central defection to plate thickness ratios 
ranging from 0.5 to 2.0 are compared with the deflected shapes obtained 
by using the plate-membrane coupling model. Next the load-displacement 
curve obtained by using the finite element model is compared with the 
load-displacement curve obtained earlier by using the plate-membrane 
coupling model. 

The transverse displacements w obtained by using the finite element 
model and the plate-membrane coupling model (Chap. 4) for various 
central deflection to plate thickness ratios ranging from 0.5 to 2.0 are 
presented in Table 5.2. The results obtained by the finite element 
model are in excellent agreement with those obtained earlier by using 
the plate-membrane coupling model. 

Figure 5.4 presents a comparison between the load-displacement 
curves obtained by using the finite element model and the plate-membrane 
coupling model. Figure 5.4 indicates that the results obtained by using 
the finite element model compare well with those obtained by using the 
plate-membrane coupling model. Figure 5.4 also shows the effect of 
nonlinearity on the plate deflection. For any given load the nonlinear 
displacement is much smaller than the linear displacement. 

As pointed out in the introduction, in the current finite element 
formulation, the membrane effects due to the mid-plane stretching are 
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2 
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Fig. 5.4 Load-Displacement Curve 



125 


considered. So it is interesting to compare the radial displacements 
obtained by using the finite element model, with those obtained by using 
the classical solution (Appendix B). Figures 5.5 and 5.6 show a com- 
parison of the radial displacements for the central deflection to plate 
thickness ratios of 0.5 to 2.0. The functional form of the radial 
displacement curve obtained by using the classical solution is identical 
for the central deflection to plate thickness ratios of 0.5 and 2.0. On 
the other hand, the functional form of the radial displacement curve 
obtained by using the finite element model and the plate-membrane 
coupling model is different for the central deflection to plate 
thickness ratios of 0.5 and 2.0. 

Figures 5.5 and 5.6 also show that the magnitudes of the radial 
displacements are three orders less as compared to the magnitudes of the 
transverse displacements. Although the radial displacements are much 
smaller in comparison with the transverse displacements, their contri- 
bution to the membrane radial strains, and the tangential strains is 
significant. These membrane strains are particularly important for 
strain energy calculations in the back face spalling model. Thus for 
accurate predictions of stresses and strains in the plate, the radial 
displacements should be incorporated into the finite element formula- 
tion. 

Figure 5.6 indicates that, for the central deflection to plate 
thickness ratio of 2.0, radial displacements obtained by using the 
classical solution are in good agreement with those obtained by using 
the finite element solution and the plate-membrane coupling model 
solution. On the other hand. Fig. 5.5 shows that, for the center 

deflection to plate thickness ratio of 0.5, radial displacements 
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Fig. 5.5 Comparison of Radial Displacements for the Center Deflection to the Plate 
Thickness Ratio Equal to 0.5 


Finite Element Solution 
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Fig. 5.6 Comparison of Radial Displacements for the Center Deflection to the Plate 
Thickness Ratio Equal to 2.0 
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obtained by using the classical solution differ considerably from the 
corresponding finite element solution and the pi ate -membrane coupling 
model solution. Therefore the functional form of the radial displace- 
ments assumed in the classical solution would predict correct radial 
displacements only for higher values of the center-deflection to plate 
thickness ratios, while for lower values of the center deflection to 
plate thickness ratios, the functional form of the radial displacements 
assumed in the classical solution would predict incorrect values of the 
radial displacements. 


5.9 Concluding Remarks 

In this chapter, circular quasi-isotropic composite plates were 
analyzed to study the large deformation behavior under point loads. A 
finite element model was formulated to study the large deformation 
behavior of these plates. A finite element model, in contrast to those 
in the literature, considers both radial and transverse displacements in 
its formulation. This model uses a linearization technique with an 
iterative procedure to obtain the large deformation shapes for thin 
composite laminates under the point loads. 

The deformation shapes obtained by using the finite element model 
agreed very well with the deformation shapes obtained by using the 
plate-membrane coupling model. The load-deflection curve was obtained 
using the finite element model and compared with the load-deflection 
curve obtained by using the pi ate -membrane coupling model. Excellent 
agreement was observed between the two results. The classical solution 
[27], based on the energy method, assumes that the functional form of 
the radial displacements curve is identical for the various values of 
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the central displacements to plate thickness ratios, whereas the finite 
element model indicates that the functional form of the radial displace- 
ments curve is different for the various values of the central displace- 
ments to plate thickness ratios. The radial displacements are found to 
be about three orders of magnitude less, as compared to the transverse 
displacements. The radial displacements obtained by using the finite 
element solution compares well with the classical solution for the 
central deflection to plate thickness ratio of 2.0, but they are 
significantly different from the classical solution for the center 
deflection to plate thickness ratio of 0.5. 
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Table 5.1 - Number of Iterations Required to Obtain Large 

Deformation Solution of a Clamped Circular Plate 
Under Central Point Load 


Plate center deflection Number of 

w o 

to thickness ratio (— iterations 


0.5 


4 


1.0 9 

1.5 17 


2 . 0 - 


29 
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Table 5.2(a) - Comparison of the Transverse Displacements w 

Obtained by Using the Finite Element Mo^el and 

the PI ate -Membrane Coupling Model for (— 2 -) Ratios 
of 0.5 and 1.0. n 


r/a 


Finite Element 
Model 


Plate-Membrane 
Coupling Model 




0.0 

0.50000 

0.50000 

0.1 

0.46810 

0.46764 

0.2 

0.40963 

0.40842 

0.3 

0.34011 

0.33791 

0.4 

0.26712 

0.26484 

0.5 

0.20418 

0.20311 

0.6 

0.13558 

0.13506 

0.7 

0.07916 

0.07762 

0.8 

0.03763 

0.03828 

0.9 

0.00983 

0.01001 

1.0 

0.00000 

0.00000 

0.0 

1.00000 

1.00000 

0.1 

0.93167 

0.93050 

0.2 

0.81098 

0.80331 

0.3 

0.67184 

0.66745 

0.4 

0.53024 

0.52491 

0.5 

0.39756 

0.39392 

0.6 

0.27246 

0.27028 

0.7 

0.16224 

0.16570 

0.8 

0.07546 

0.07719 

0.9 

0.02024 

0.02074 

1.0 

0.00000 

0.00000 
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Table 5.2(b) - Comparison of the Transverse Displacements w 

Obtained by Using the Finite Element Mo^jel and 

the PI ate -Membrane Coupling Model for (-£■) Ratios 
of 1.5 and 2.0. n 



r/a 

Finite Element 

Plate-Membrane 



Model 

Coupling Model 

( nr ) - ‘- 5 

0.0 

1.50000 

1.50000 

0.1 

1.38780 

1.38595 


0.2 

1.19976 

1.19548 


0.3 

0.99214 

0.98512 


0.4 

0.79318 

0.78421 


0.5 

0.59440 

0.58841 


0.6 

0.40648 

0.40245 


0.7 

0.24651 

0.25309 


0.8 

0.11811 

0.12140 


0.9 

0.03186 

0.03286 


1.0 

0.00000 

0.00000 

=£>* 

w 

ro 

o 

0.0 

2.00000 

2.00000 

0.1 

1.83683 

1.83399 


0.2 

1.57758 

1.57111 


0.3 

1.30415 

1.29297 


0.4 

1.03569 

1.02026 


0.5 

0.78009 

0.77129 


0.6 

0.54370 

0.53813 


0.7 

0.33418 

0.34512 


0.8 

0.16317 

0.16870 


0.9 

0.04395 

0.04561 


1.0 

0.00000 

0.00000 




Chapter 6 


PLATE MECHANICS EXPERIMENTS 

6.1 Introduction 

In Chaps. 4 and 5, the plate-membrane coupling model and the finite 
element model were developed to study the large deformation behavior of 
circular quasi-isotropic plates. Both of these models assume that the 
quasi-isotropic T300/5208 graphite/epoxy laminates with stacking 
sequence [45/0/-45/90] s , are flexurally isotropic. In this chapter this 
assumption of flexural isotropy is experimentally verified. 

The large deformation shapes and the load-displacement curves for 
the circular quasi-isotropic plates under central point loads were 
obtained by conducting series of tests on the circular quasi-isotropic 
laminates. First, the experimentally obtained deflected shapes are 
compared with the analytical deflected shapes obtained by using the 
plate-membrane coupling model and the finite element model. Then, the 
experimental load-displacement curve is compared with those of the 
plate-membrane coupling model and the finite element model. 

6.2 Test Set-up 

To study the large deformation behavior of circular plates under 
central point loads, a torque bolt arrangement and a servo-hydraulic 
testing machine were used. The torque bolt arrangement was used to 
obtain the large deformation shapes of the circular plates under central 
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point loads and the servo-hydraulic testing machine was used to obtain 
the load-displacement behavior of the circular plates under central 
point loads. For both tests, quasi-isotropic graphite/epoxy laminates, 
with the stacking sequence of [45/0/-45/90] s and the material properties 
given in Table 2.1 were used. Circular laminates with 76.2 mm (3 
inches) or 101.6 mm (4 inches) diameters were bolted between two steel 
annular plates as shown in Fig. 6.1. These steel plates provided 

clamped boundary conditions for the specimens. A torque bolt 

arrangement was used to push a 25.4 mm diameter steel ball against the 
laminate, to provide the desired static loading (Fig. 6.2). Although 

the ball diameter was 25.4 mm, the contact radius between the ball and 
the plate was very small (of the order of 1 mm) and thus pushing the 
steel ball simulated the desired static point load conditions. 

A traversing horizontal direct current differential transducer 
( DCDT) and a vertical DCDT were used to measure the deflected shapes of 
the clamped circular quasi-isotropic laminates as shown in Fig. 6.1. 

To obtain a load-deflection curve, the static loading tests were 
conducted in a servo-hydraul ic testing machine. The test specimens were 
clamped to a platform, shown in Fig. 6.3, which was mounted on the 

hydraulic ram, and load was applied to the center of the specimen by 
means of a punch, tipped with a 25.4 mm diameter steel ball. A DCDT 
displacement gauge was attached to the indenter to measure plate deflec- 
tion as the load was applied. The electric signals from the load cell 
and the displacement gauge were given as an input to the x-y plotter to 
allow direct load-displacement plotting. 
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Hydraulic Ram 


g. 6.3 Static Loading on the Quasi-Isotropic Circular 
Laminate 
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6.3 Test Data 

Circular quasi-isotropic plates of two sizes With diameters 76.2 mm 

and 101.6 mm were tested. By using the torque bolt arrangement, the 

76.2 mm diameter plate was loaded till the desired center deflection was 

obtained. This center displacement was measured by using the vertical 

DCDT. Once the desired center displacement was obtained, the first set 

of observations were made along the fiber direction (0 = 0°) of the 

eighth ply (farthermost ply from the point of load application). The 

second and third set of observations were made on a line perpendicular 

to the fiber direction of the eighth ply, (0 = 90°) and on a line 45° 

(0 = 45°) to the fiber direction of the eighth ply, respectively. 

These three sets of observations were recorded for central deflection- 

w o 

to-plate-thickness ratio (— p-) of 0.5, 1.0, 1.5, and 2.0. Similar 

observations were made for 101.4 mm diameter plate. 

The load-displacement data were obtained by using the servo- 
hydraulic machine. The 76.2 mm diameter and the 101.4 mm diameter 
plates were loaded and the corresponding central deflections were 
recorded till the central deflection was about two times the plate 
thickness. 


6.4 Results and Discussion 

In this section, first the validity of flexural isotropy in 
circular quasi-isotropic laminate is discussed. Then the measured 
deflected shapes are compared with the plate-membrane coupling model and 
the finite element model. 

To verify the flexural isotropy in the circular quasi -i sotropic 
laminates, experimental measurements made along three diametral lines at 



139 


6=0°, 0 = 45° and 6 = 90°, for various values of the central 

deflection-to-plate thickness ratios ranging from 0.5 to 2.0 were 

w 0 

compared. Figure 6.4 presents the typical deflected shapes for (— ^-) 
ratios of 0.5 and 1.5 along 0 = 0° and 90° radial lines. The 
deflected shapes for ratios of 1.0 and 2.0 showed similar trends and 
hence are not shown. The experimental observations showed that the 
deflected shapes along 0 = 45° were bounded by the shapes along 
0=0° and 0 = 90°. The maximum differences between the 0 0° and 
90° observations occured near (-£•) = 0.4 and were about six percent of 

a 

the maximum center deflections. Thus the experimental observations 
indicate that clamped circular quasi -i sotropic laminates, subjected to 
axisymmetric loading, exhibit near flexural isotropy and show nearly 
axisymmetric bending behavior. 

Figures 6.5 and 6.6 present the measured deflected shapes at 

0=0° and the analytical deflected shapes obtained with the plate- 

w 

membrane coupling model for ratios ranging from 0.5 to 2.0. 

Since the deflected shapes obtained by using the finite element model 
are almost identical to those obtained by using the plate-membrane 
coupling model (Table 4.1), only the plate-membrane coupling model 
results are compared with the experimental data. In Figs. 6.5 and 6.6 
the experimental results are shown as circular symbols and solid lines 
represent the results of the pi ate -membrane coupling model. Figures 6.5 
and 6.6 show excellent agreement between analytical and experimental 
results. 

Figure 6.7 presents the measured load-deflection curve and the 
analytical load-deflection curve obtained by using the plate-membrane 
coupling model and the finite element model. The analytical results 


t 



Fig. 6.4 Axisymmetric Validation for Quasi-Isotropic 8-Ply T300/5208 Laminate 




142 



Fig. 6.6 Deflected Shapes for the Laminate Under Quasi- 

Static Point Loads for w /h = 1.5 and 2.0 
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were found to be within six percent of the experimental results. These 
differences between the analytical and experimental results are probably 
due to experimental errors which might have been introduced while 
practically simulating the clamped boundary conditions. 

6.5 Concluding Remarks 

Eight ply T300/5208 quasi-isotropic circular plates were tested to 
verify the assumption of flexural isotropy and to study the large 
deformation behavior under central point loads. The experimental 
observations indicate that clamped circular quasi-isotropic laminates, 
subjected to a central point loading, behave as if they are almost 
flexurally isotropic and they show nearly axisymmetric bending behavior. 

The experimentally obtained large deformation shapes and the load- 
displacement curves were compared with those obtained with the plate- 
membrane coupling model and the finite element solution. The experi- 
mental results are in good agreement with the analytical results 
obtained by using the plate-membrane coupling model and the finite 
element model . 

In the next chapter, the large deformation shapes of the quasi- 
isotropic circular plate obtained by using the plate-membrane coupling 
model and the finite element model are used in conjunction with the 
fracture mechanics approach to describe the back face spalling 
phenomenon in thin quasi-isotropic circular laminates. 



Chapter 7 


BACK FACE SPALLING MODEL 
7.1 Introduction 

In Chaps. 4 and 5, the plate-membrane coupling model and the finite 
element model were developed to study the large deformation behavior of 
clamped circular quasi-isotropic laminates under static point loads. 
These models were verified by conducting tests on the clamped circular 
laminates. In this chapter a simple analytical model is developed to 
describe the back face spalling phenomenon in thin quasi-isotropic 
laminates. The model is based on the fracture mechanics principles and 
uses the large deformation shapes of the plates obtained earlier by 
using the plate-membrane coupling model. This analytical model is 
verified by conducitng the experiments on thin quasi-isotropic T300/5208 
circular plates. 


7.2 Back Face Spalling Mechanism 

As discussed in Chap. 1, for quasi-isotropic laminates under the 
central point loads, first visible damage occurs in the farthest ply 
from the load (i.e. 8th ply). The damage is in the form of a matrix 
cracking parallel to the fibers. Further increase in load level results 
in two dominant cracks in the 8th ply and surrounding these two cracks 
delaminations form between 7th and 8th plies. Additional increase in 
load causes these delaminations and cracking to propagate ustably in the 
fiber direction of the bottom ply. This unstable growth of delamina- 
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tions and the propagation of the two dominant cracks is known as the 

back face spalling phenomenon. During the back face spalling action, 
curved trajectory of the 8th ply center strip (formed by the two 
dominant cracks in the 8th ply) running under the load point in the 
fiber direction changes virtually to a straight line between the load 
point and the support (Fig. 7.1). A large fraction of strain energy 

from this 8th ply center strip is shed during this spalling action. 

This mechanism of the back face spalling is modeled by using the 
concepts of fracture mechanics in the following sections. 

7.3 Fracture Mechanics Concepts 

Any general deformation of a cracked body can be described by 
combination of three independent modes of deformation (Fig. 7.2). The 
three independent modes are: opening mode, sliding mode and tearing 

mode. Normal stresses give rise to the opening mode denoted as mode 

I. Inplane shear results in mode II or sliding mode. The tearing mode 
or mode III is caused by out-of-plane shear. 

7.3.1 The Griffith Criterion 

One of the basic criterion for fracture was established in 1921 by 
Griffith [54]. The Griffith energy criterion for fracture can be stated 
as "crack growth can occur if the energy required to form additional 
crack surfaces can be delivered by the system." 

For example, consider a cracked plate of thickness B under a load 
P, as shown in Fig. 7.3. For a crack of size a the linear load- 
displacement relationship can be represented by line OA in Fig. 7.4. 
Similarly for a crack of size a+da the load-displacement relation can 





Fig. 7.1 Back Face Spalling in a Composite Plate 
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Fig. 7.2 Three Modes of Fracture 
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be represented by the line 0E (the line OE has a lower slope than the 
line OA because of reduced stiffness of the plate). If the ends of the 
plate are fixed and if crack extension takes place from a to a+da at a 
load Pj, the displacement remains constant and the load drops from point 
A to point B (Fig. 7.4). During this process elastic energy, repre- 
sented by the area of triangle OAB, is released, which is A P.Y. 

If crack extension takes place at constant load the displacement 
increases from point A to point E, i.e. by an amount AY (Fig. 7.4). 
Therefore the work done W, by the load is 

U = P 2 • AV (7.1) 

The increase in the elastic energy of the plate is 

area (OEF) - area (OAC) = -^ P^v+Av) - ^ P^*V = -^ P^-AV (7.2) 

This increase in energy has to be provided by the load. The energy 
provided by the load is area (AEFC) 

area (AEFC) = P L • AV (7.3) 

Equations (7.2) and (7.3) indicate that there remains an amount of 
energy equal to area (OAE). Area (OAE) can be expressed as 

area (OAE) = ^ P^AV = AP*V + AP*AV 


(7.4) 
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Fig. 7 . 4 


Load-Displacement Curve 
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In Eq. (7.4) -Jj- AP*AV is the second order term and hence can be 
neglected. Thus it follows that area' (OAB) => (OAE). This means that 
energy available for crack growth is the same in fixed grip and constant 
load case. 

As seen earlier in the case of fixed grips the elastic energy is 
released as the crack grows from a to a+da. This elastic energy 
released per unit crack extension is called "elastic energy release 
rate" and is denoted by G. Since elastic energy is proportional to the 
strain energy, G is usually expressed in terms of the strain energy and 
is called "strain energy release rate." 

7.3.2 Critical Strain Energy Release Rate G 

Griffith [54] postulated that the crack growth occurs only when the 
strain energy release rate G exceeds certain critical values. This 
critical value is called a "critical strain energy release rate" and is 
denoted by G c . This quantity G c is a material property like the yield 
stress. Young's modulus and is different from material to material. 

7.4 Back Face Spalling Model 

As discussed earlier in the section 7.2, in quasi-isotropic 
laminates under the point loads, first visible damage occurs in the 
lowest ply in the form of two dominant cracks and surrounding these two 
cracks delaminations form between the 7th and 8th plies. With further 
increase in load at a certain plate center deflection w Q , massive back 
face spalling occurs. During this spalling action the 8th ply center 
strip (formed by the two dominant cracks in the 8th ply) which is 
previously deformed in the same shape as that of a circular plate 
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changes virtually to a straight line between the load point and the 
support. This failure mechanism suggests that the complete back face 

spalling event can be explained by modeling the center strip in the 8th 

ply. 

The center strip model developed here to describe the back face 

spalling phenomenon is based on the following assumptions: 

1. Before the back face spalling occurs, the center strip in the 

8th ply has the same deformed shape as that of a quasi-isotropic 
plate under the point load. Thus the radial displacements (u), 
the transverse displacements (w), and the rotations (^r)> in 
the center strip can be obtained by performing the large deflec- 
tion analysis of the quasi-isotropic circular plate under point 
loads. 

2. The initial delamination between the 7th ply and the center 

strip in the 8th ply is assumed to exist at the radius of 

inflection (the radius at which curvature in the strip is zero). 

3. The change in stiffness of the plate due to the back face 

spalling action is assumed negligible and hence is not 
considered in the center strip model. 

4. The effects of internal damage in the circular plate on the 

deformed shape of the plate are assumed to be negligible. 

5. The center displacement of the plate is held constant and it 

does not change during the back face spalling action. 


With these assumptions the center strip model was formulated and 
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solved to compute the strain energy release rates during the back face 
spalling action. 

Consider a center strip in the 8th ply. The deformed shape of the 
strip is the same as that of the circular plate. The strip has a width 
b and the thickness t (t = one ply thickness). This center strip was 
discretized into m number of elements with (m+1) nodes. The radial 
displacement (u), the transverse displacement (w), the rotation , 

and the curvature at each node were calculated by using the plate- 

membrane coupling anlysis as described in Chap. 4. An initial 

delamination was assumed to exist between the center strip (of the 8th 
ply) and the 7th ply at the radius of inflection as shown in Fig. 7.5. 
This initial delamination was assumed to have a length equal to BD as 
shown in Fig. 7.5 where nodes B and D are the nodes adjucent to the node 
at the radius of inflection. Due to the presence of this initial 
delamination, the original deformed shape of the center strip ABCDEF as 
shown in Fig. 7.5 changed to the new deformed shape ABC'DEF. The 

center strip initially deformed between the nodes BD as BCD became 
straight as BC'D. Due to this change the original displacements, 

rotation and curvature at the node C have changed. The new displace- 
ments, rotation and curvature at the node C were calculated as follows: 
The radial displacement at the node C' equal to 

( u n ' u R ) 

U C = ~2 (Ar r ' (Ar) + U B <7 ' 5) 


The transverse displacement at the node C' equal to 


NODE AT RADIUS OF AT Ar AT 



Fig. 7.5 Delammations in Center Strip 
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w 


C 


(Ar) + Wg 


(7.6) 


The rotation at the node C' equal to 


/-dW’i _ ^ W D 

l 3F J C' T - 


W p ) 

TaFT 


(7.7) 


The curvature at the node C‘ using finite difference approximation at C* 
is 


r d 2 w A _ W D " 2 V + W B 

l 7 ) . 2" • ('Ar) 


(7.8) 


dr C‘ 


Substituting equation (7.6) for Wq* into Eq. (7.8) 


(^) - 0 


(7.9) 


dr C* 


Thus with an initial delamination the center strip had the same 
nodal displacements, rotations, and the curvatures as that of a plate, 
except at the node C' where the displacements and rotation are given by 
Eqs. (7.5)— (7.8). The curvature at the node C', however, is zero. By 
using these values of displacements, rotations and curvatures strain 
energy of the center strip was calculated as follows: 

The strain energy U of a body in terms of stress (a) and strain 
( e) is given by: 



a e d(vol ) 


(7.10) 
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Only the transverse displacements (w) and the radial (axial) displace- 
ments (u) were assumed to be present in the strip. With - these 
assumptions, the strip was modeled as a beam and hence only one stress 
o r exists in the center strip. The strain energy of the strip 

idealized as a beam is, 

U d (vol) (7.11) 

c vol r 


where vol is the volume of the center strip and E rr = E^ is the 
elastic modulus of the strip in the fiber direction. 

In Eq. (7.11) e r is the radial (axial) strain and is the sum of 
the membrane strains and the bending strains: 


e 

r 


du 1 fdw^2 

' ¥ 7 l 3r j 



(7.12) 


where z in the Eq. (7.12) was measured with respect to the plate 
coordinate system. 

A substitution of Eq. (7.12) in Eq. (7.11), gives the total strain 
energy of the strip U as 

E 2 2 

ii 11 ( r d cj 1 /-dw-i2 d w . » in n 

U = — * — J U- + -t UrrJ " 2 7 ) d (vol) (7.13) 

Since z was measured with respect to the plate coodinate system, 
the limits of integration for the 8th ply center strip in z direction 
are (3h/8) to (h/2), therefore 
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U = 


E ll b 


J* / 


h/2 


3h/8 


f du + 1 

W 7 


/dw^ 

W J 


d 2 w^2 . . 

z — _) dr dz 

dr 


(7.14) 


In order to evaluate the integral in Eq. (7.14), the center strip 
was discretized into m number of elements. The nodal displacements, 
the roation, and the curvatures were computed by using the plate- 
membrane coupling model. These displacements, rotation and curvature 
were modified at the node C* (Fig. 7.5) by using Eqs. (7.5)-(7.9). 
Using these modified nodal displacements, Eq. (7.14) was numerically 
integrated to obtain a strain energy U of the center strip of the 8th 
ply. This procedure was repeated for various values of m, where number 
of elements m, were chosen as multiples of 10. About 60 elements and 90 
elements were found necessary for a 25.4 mm and 38.1 mm length center 
strip respectively, to yield a converged solution for the strain energy 
U. Once the strain energy U was computed a delamination of one element 
size was allowed to grow in the outward direction as shown in Fig. 7.5. 
Due to this growth of delamination the deformed shape of the strip 
ABC'DEF changed into a new deformed shape ABC'‘D'EF as shown in Fig. 
7.5. Thus due to growth of delamination the deformed shape of the strip 
between the nodes B and E is a continuous straight line. The new 
displacements, rotation and the curvatures at nodes C 1 ', D' were calcu- 
lated as follows: 

The radial displacements of the nodes C ’ * and O’ are 


(u E ~ u B ) 

U C" ° -3-T • “■ + U B 


(u, - u R ) 

V = TW * 2(Ar) + U B 


(7.15) 


(7.16) 
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The transverse displacements at the nodes C' ' and O' are 

(W — W rj ) 

V g 3 • (ArT • Ur) * w 8 (7 - 17 > 

“O' ' T^TZFT ‘ 2 f4r) + “B (7 - 18) 


The rotations at the nodes C‘‘ and O’ are 


r dw >, _ cdw^ _ (w E " W B* 

^dr'C' 1 ~ '■(JF-'D 1 3 • (A r ) 


(7.19) 


Since the deformed shape of the strip between the nodes B and E was 
straight, the curvatures at the nodes C" and O' were zero. Therefore 

= (^) =0 (7.20) 

dr C' 1 dr O' 

Thus for the new deformed shape ABC‘'D'CF except at the nodes C’ 1 
and D' the nodal displacements, rotations, and curvatures of the center 
strip are the same as those obtained by using the plate-membrane 
coupling model. For the nodes C 1 ', D', the displacements, rotations, 
and curvatures are given by Eqs. (7. 15) -( 7.20) . Using these values of 
nodal displacements, rotations, and curvatures into equation (7.14), the 
new strain energy U1 of the center strip with the deformation shape 
ABC'D’EF was calculated. 

This strain energy U1 corresponding to the deformed shape ABC'D'EF 
was found less than the strain energy U corresponding to the deformed 
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shape ABC'DEF. Thus during the delamination growth from BD to BE (Fig. 
7.5), the amount of strain energy released All - is given by 


A U = U - U1 


(7.21) 


Expressing the strain energy released All in terms of strain 
energy release rate G, one has 


G = 


au AU 

bJ " A* 


(7.22) 


where 


A A = b x Ad Q (7.23) 

b = width of the strip 
Ad Q = length of the delamination growth 
This G corresponds to a delamination length of d Q + Ad Q . 

By using the above procedure, the strain energy release rates were 

calculated by incrementing the delamination outward by one element 

length each time, until the delamination reached the clamped edge (node 

F, Fig. 7.5). This strain energy release rate analysis was repeated for 
w 

several (-^) ratios and for two plate sizes with radii 25.4 mm and 
38.1 mm. 


7.5 Strain Energy Release Rate Results 

In this section the strain energy release rates obtained by using 
the center strip model are presented for two plate sizes with radii 25.4 
mm and 38.1 mm. For each plate size strain energy release rates were 
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determined for each delamination length. Figures 7.6 and 7.7 present G 
as a function of delamination length for various values of plate center 
deflection w Q . For all delamination lengths. Fig. 7.6 shows that larger 
values of center deflection w Q gave larger values of the strain energy 
release rates. Figure 7.6 also indicates that beyond delamination 
length of 15 mm, the strain energy release rates reach a constant value 
and do no depend on the delamination length. 

Figure 7.7 presents the strain energy release rate variation with 
delamination length for a plate radius of 38.1 mm, for various values of 
center deflection. This figure shows similar trends as Fig. 7.6. 

7.6 Discussion of the Back Face Spalling Phenomenon 

As pointed out earlier, the back face spalling pehnomenon occurs 
spontaneously (i.e. a center strip from the 8th ply peels off unstably). 
This unstable peeling action suggests that, during delamination growth 
between the 7th ply and the center strip from the 8th ply, the strain 
energy release rate G must exceed the critical strain energy release 
rate for delamination growth. Assuming that the spalling occurs due to 
the peeling action alone, i.e. only mode I is present, then by using the 
cirical strain energy release rate value of the T300/5208 graphite/epoxy 
material in the opening mode (Gj) c , the back face spalling damage can be 
predicted by using the center strip model. 

Consider a 25.4 mm radius plate. The critical strain energy 
release rate value (Gj) c , for the T300/5208 graphite/epoxy material is 
0.1 KJ/m^ [55]. Figure 7.6 shows that for the plate center deflections 
of 1.60 mm or less, the calculated G values to grow initial delamination 
length BC'D to new delamination length BC'D'E (Fig. 7.5), are less than 
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0.1 KJ/rrr or {G j ) c value. This indicates that for the plate center 

deflections of 1.60 mm or less, the initial delamination between the 

center strip and the 7th ply will not grow. However, for the plate 
center deflections of 1.70 mm or more. Fig. 7.6 shows that the 

calculated G values to grow the initial delamination length BC'D to new 
delamination length BC'D'E (Fig. 7.5), exceed ( G j ) c value. Therefore, 
for plate center deflections equal to or greater than 1.7 mm the initial 
delamination between the center strip and the 7th ply will grow unstably 
and the center strip will peel off, thus resulting into the back face 
spalling action. Figure 7.6 also indicates that if the initial 

delamination lengths are large, the back face spalling action will occur 
at lower values of central deflection w Q . For example consider an 

intitial delamination length of 5 mm. From Fig. 7.6 the calculated G 
value required for growth in the initial delamination length of 5 mm 
exceeds 0.1 KJ/m^ or ( G j ) c value when w Q is 1.60 mm. That means even at 
the center deflection of 1.60 mm, the initial delamination of 5 mm 
length will grow unstably resulting in the back face spalling action. 

Next consider the larger plate of 38.1 mm radius. Figure 7.7 shows 
that for the plate center deflection of 2.6 mm or more, the calculated G 
values to grow an inital delamination length BC'D to new delamination 
length BC'D'E (Fig. 7.5), esceeds (Gj) c value. Therefore, for plate 
center deflections equal to or greater than 2.6 mm the initial delamina- 
tion between the center strip and the 7th ply will grow unstably and the 
center strip will peel off, thus resulting into the back face spalling 
action. Figure 7.7 also indicates that if the delamination lengths are 
large, the back face spalling action will occur at lower values of 

o* 


center deflections w. 



7.7 Back Face Spalling Experiments 


To verify the back face spalling model, experiments were performed 
on clamped circular quasi-isotropic plates. Quasi-isotropic 

graphite/epoxy laminates, with the stacking sequence of [45/0/-45/90] s 
and the material properties given in Table 2.1 were used for the tests. 
Circular laminates with 25.4 mm and 38.1 mm radii were bolted between 
two annular steel plates. By using the torque-bolt arrangement (similar 
to the one described earlier in Chap. 6), a 25.4 mm diameter steel ball 
was pushed against the laminate. Although the ball diameter was 25.4 
mm, which was half the size of the plate diameter, the contact radius 
between the ball and the plate was very small (of the order of 1 mm) and 
thus pushing of the steel ball simulated the desired static point load 
conditions. A vertical DCDT was used to measure the center deflection 
w Q . The static load on the plate was gradually increased by using the 
torque bolt, till the back face spalling occurred. At this instant the 
center defleciton w Q of the plate was measured. From the experiments it 
was observed that for 25.4 mm radius plate, the back face spalling 
occurred within the center deflections, ranging from 1.80 mm to 1.90 mm. 
For 38.1 mm radius plate the back face spalling occurred at higher 
center deflections, with the values ranging form 3.00 mm to 3.1 mm. 

The back face spalling experiments show that for the two plate 
sizes studied, a plate with 25.4 mm radius failed at a mean center 
deflection of 1.85 mm. The critical strain energy release rate (G c ), 
corresponding to the center deflection w Q = 1.85 mm from Fig. 7.6 is 

O 

about 0.14 KJ/m . In the case of 38.1 mm radius plate back face 
spalling occured at a mean center deflection of 3.05 mm. The critical 
stran energy release rate G corresponding to the center deflection w 

o 
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= 3.05 mm form Fig. 7. 7 is 0.147 KJ/m^. Thus the experimentally obtained 
critical strian energy release rates G c , are about 0.14-0.15 KJ/m^. 
Comparing these values with the opening mode critical strain energy 
release rate (Gj) c [55], G c is somewhat higher than { G j ) c - This may be 
due to the presence of shear sliding mode (mode II) in addition to the 
opening mode (mode I) during the back face spalling action. 

Earlier investigators [56,57] have shown that for brittle systems 
like T300/5208 graphite/epoxy material under the mixed mode situations 
(opening mode and shear sliding mode), the failure is predominantly 
governed by the critical strain energy release rate in the opening mode. 
The present experimental values of G c are very close to (Gj) c . Hence, 
even though in the back face spalling phenomenon a mixed mode (mode I + 
mode II) situation may exist, the critical strain energy release rate in 
the opening mode (Gj) c , probably governs the back face spallng mechanism 
in thin composite laminates. 

7.8 Prevention of Back Face Spalling in Thin Composite Laminates 

The analytical and experimental results form the back face spalling 
model suggest that, the composite laminates with higher (Gj) c values 
than T300/5208, will have better resistance to the back face spalling. 
For example consider composites AS1/3501-6 made by Hercules. Russel and 
Street [58] obtained the (Gj) c value of AS1/3501-6 as 0.15 KJ/m^. This 
value is fifty percent higher than that for T300/5208 composites. 
Therefore AS1/3501-6 composites may have better resistance to the back 
face spallng when compared to T300/5208 composites. Recently consider- 
able attention is focused on the development of tough resins. The tough 
systems of composites usually have higher (Gj) c values than brittle 
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system of composites [56]. Mall and Johnson [59] determined the ( G j ) c 
values for FM-300 tough matrix composite by Hexel. They found that 
(Gj) c values are 7 to 10 times higher than the brittle systems (i.e. 
T300/5208 and AS 1/3501—6) - Hunston [60] and O'Brien et al. [61] 
obtained (Gj) c values for F-185 matrix composite by American Cynamide. 
They reported the value of (Gj) c for T6000/F185 composite as 2 KJ/m^, 
which is about 20 times higher than the T300/5208 values. Thus tough 
systems of composites may have superior resistance to back face spalling 
when compared to brittle systems of composite. 

Some investigators [62] have found that "stitching" of the laminate 
decreases delamination. The stitching of the laminate may be 

particularly useful to prevent the back face spalling, as stitching has 
a restraining effect on peeling action. Another possible way to prevent 
the back face spalling, is to use woven ply on the back of the laminate, 
because the weave will not allow center strip from the bottom ply to 
peel off without involving cross-ply strips. 

Finally, in the present analysis the plates were assumed to be 
stress free. If the plates are pre-stressed, the plates would react 

differently to the center point load. For the center point load, the 
center deflection w Q of the tensile pre-stressed plate would be less 
compared to the center deflecton w 0 of the plate with zero pre- 
stresses. This is because the tensile pre-stressed plate is more stiff 
compared to the plate with zero pre-stresses [27]. This suggests that 
the back face spalling would occur at higher loads in thin composite 
plates pre-stressed with tensile stresses compared to the initially 
stress free plates. 
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7.9 Conclusions 

In this chapter a simple analytical model was developed to 
understand the back face spalling phenomenon in thin quasi-isotropic 
laminates. The model uses the large deformation plate shapes obtained 
by using the plate-membrane coupling model in conjunction with the 
fracture mechanics principles. The strain energy release rates during 
back face spalling action were obtained by using the analytical model 
for two plate sizes with plate radii of 25.4 mm and 38.1 mm. Experi- 
ments were performed on the 8 ply T300/5208 quasi-isotorpic circular 
plates to obtain critical strain energy release rates during back face 
spalling action. Good agreement was observed between experimental and 
analytical results. Experimental and analytical results indicated that 
the back face spalling phenomenon in 8 ply T300/5208 quasi-isotropic 
laminates is governed by the critical strain energy release rate ( G j ) c 
in the opening mode. Therefore the back face spalling in composites can 
be reduced by using the material with higher (Gj) c values than ( G j ) c 
value of the brittle systems composite like T300/5208. Another possible 
way to prevent the back face spalling, is to stitch the laminate or to 
use a woven ply on the back of the laminate. 
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Chapter 8 

SUMMARY AND CONCLUSIONS 

The present investigation is concerned with the failure analysis of 
composite plates under low-velocity impact loads. This analysis is 
particularly important to understnad the damage mechanics and to develop 
composites with improved impact resistance. The sate-of-art of the low- 
velocity impact problem was reviewed and important aspects of the 
problem such as deformation mechanics, back face spalling phenomenon 
were investigated. Three models, plate-membrane coupling model, finite 
element model, and back face spalling model were developed. These 
investigations should help in understanding the low velocity impact 
problem and lead to better design and analysis of composite structures. 

The low velocity .impact problem in T300/5208 graphite/epoxy 
circular quasi-isotropic plates was formulated by replacing impact type 
point loads with equivalent quasi-static loads. The quasi-isotropic 
plates were modeled as isotropic plates having the flexural stiffness 
components equivalent to the flexural stiffness components of the quasi- 
isotropic plates. A plate-membrane coupling model and a finite element 
model were developed to obtain the large deformation behavior of 
circular composite plates under point loads. These models were verified 
by conductng plate mechancis experiments. From the study of these 
analytical models and experiments the following conclusions can be made: 
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1. A circular quasi-isotorpic plate with stacking sequence of 
[45/0/-45/90] s -can be modeled as an isotropic plate having 
flexural stiffness components equivalent to those of a quasi- 
isotropic plate. 

2. The large deformation behavior of circular isotropic membranes 
subjected to arbitrary axisymmetric loading can be obtained by 
solving a single nonlinear governing equation in terms of radial 
stress using the finite difference method in conjunction with 
Newton-Raphson method. 

3. Accurate large deformation behavior of thin circular quasi- 
isotropic plates can be obtained by using a simple plate- 
membrane coupling model. 

4. The functional form of the deformed shape of the plate 
undergoing large deformations is different from the small 
deflection plate solution. The deformed shape is a function of 
the center point displacements and thus is different for 
different load levels. Furthermore, for a plate undergoing 
large deformations, there is an inward movement of the radii of 
points of inflection. 

In addition to the deformation mechanics, a study was undertaken to 
develop a damage mechanics. From the experiments on thin laminates it 
was observed that first visible damage in circular quasi-isotrpic 
laminates under quasi-static point loads occur on the back surface of 
the laminate in the form of spalling. To understand this back face 
spalling phenomenon, a simple model using the large deformation behavior 
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of thin composite plates in conjunction with fracture mechanics 
principles was developed. This back face spalling model was verified by 
conducting experiments on thin circular quasi -i sotropic plates. The 
study resulted in the following conclusions: 

1. The back face spalling action in thin composite laminates is a 
spontaneous action and can be modeled by using the fracture 
mechanics principles. 

2. Mixed mode (I + II) type deformations probably occurs during the 
back face spalling action. However, analysis and experiments 
suggest that the mode I (opening mode) may be the dominant 
mechanism and governs the delamination growth during the 
spalling action. 

3. Back face spalling model shows that the back face spalling 
occurs only when the strain energy release rate G, exceeds the 
critical value of strain energy release rate in the opening mode 
(Gj) c . This implies that back face spalling can be reduced or 
prevented by: 

(a) use of tough composites having higher values of (Gj) c 
compared to brittle resin systems like T300/5208. 

(b) Stitching of the laminate, as stitching has a 
restraining effect on peeling action. 

(c) To use woven ply on the back of the laminate, because 
the weave will not allow center strip from the bottom 
ply to peel off without involving cross-ply strips. 
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The present investigation has also shown that future research is 
warranted in the following areas: 

1. To include effects of the internal damages in composite plates 
and develop more general deformation mechanics. 

2. To extend the back face spalling model for rectangular composite 
plates by developing 2-D analysis. 

3. To devleop a technique for identifying individual strain energy 
release rates in different modes I, II, and III during back face 
spalling action. 

4. To determine the effects of stitching on the back face spalling 
action in thin composite laminates. 

5. To extend present analysis to thick composite plates by 
incorporating the effects due to shear. 

6. To perform a dynamic analysis and compare it with the quasi- 
static analysis. 
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APPENDIX A 

NEWTON-RAPHSON METHOD FOR NONLINEAR SYSTEM 

The purpose of this appendix is to present the Newton-Raphson 
method for solution of nonlinear simultaneous algebraic equations. 
First a method is explained with reference to a function involving one 
variable. Next, the method for a general n variable system is 

presented. Last, the method is illustrated with the set of nonlinear 
equations involving 5 variables derived in Chap. 3 for the membrane 
problem. 


A.l Newton-Raphson Method for One Dependent Variable 

Consider a function F(x) of one dependent variable x. The 
objective is to find the root of the equations F ( x ) = 0. Figure A-l 
presents the function F(x) graphically. The Newton-Raphson method is an 
iterative method which continuously updates an initial approximation 
until the actual root is found. Consider an approximation x = x Q as the 
root. The value of F(x) is F(x 0 ) at point P as shown in Fig. A-l. At P 
draw a tangent to the curve. The tangent intersects the X axis at T. 
Therefore the next approximation for the root is x^ = x Q + Ax q where 
Ax q is MT in Fig. A-l. Next draw a tangent at point Q. This tangent 
intersects the X axis at Tj_. Therefore the next approximation for the 
root is x 2 = X 1 + Ax l * Imagine a third tangent is drawn at R, this 
tangent will cut the X axis at some point T£ between Tj and S. 
Therefore points T, Tj, I 2 , ... will approach the point S as a limit. 
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that is, the intercepts OT, OT^, OT 2 , ... will approach the intercept OS 
as a limit. But OS represents the real root of the equation F(x) = 0. 
Hence the quantities OT, OTj, OT 2 , ... are successive approximations to 
the desired root. The fundamental formula for finding the root of the 
equation F(x) = 0 can be derived by using Fig. A-l as follows: 


Consider Fig. A-l. Let MT = Ax Q , and TTj = Ax^ , etc. The 

slope of the graph at P is F'(x n ). From the Fig. A-l, PM = F(x n ) and 

F (x ) 

slope at the point P = tan < X TP = F* ( x Q } = - —r- Therefore 


Ax 


- F <*„> 
r-ry- 


(A-I) 


The improved value of the root is then 


x 


1 


+ Ax 

0 


Similarly succeeding approximations are 


(A- 2 ) 


where 


where F'(x^) 


x^ = + Axj 

X3 = x 2 + Ax 2 
x n = x n-l + Ax n-1 

- F (x.) 

= (x"7)' 

is the derivative of the F{x) at . 


(A— 3) 


(A-4) 


In the one-variable 


O 3 
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case F'(x) = [ J( ) 3i x i » where [J] is called as Jacobian. 

The value of x is updated till the (a | < e, where e is the 
specified tolerance. The corresponding root of the equation F(x) = 0 is 
x.j . The above procedure can be generalized for n variables as follows: 

A. 2 Newton-Raphson Method for n Variables 

Consider a system of simultaneous nonlinear algebraic equations 
with n variables, (x lf X 2 « * 3 ,... x n ) as 


f 2 ( Xj » X 2 , x^ 


xn) = R. 


f 2 ( x 2 * x 2* x 3 •••••••• x „) R. 


‘n ’'2 


(A-5) 


f n < x l> x 2 ’ x 3 x n ) * R n 


Equations in (A-5) can be represented as 


(F (x) } = { (f 2 - R2). (f 2 ' R 2 } (f n ~ R n )} (A-6) 


The approxmate solution of the nonlinear system (Eq. (A-5)) can be 
obtained starting with an intial approximation X Q as follows: 

With the initial approximation X Q , [J (X Q )] is first evaluated 
as: 
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EJ(X 0 >] = 


af : (x) 
5x7 


5f 2 (x) 
5x7 


8f n 1x1 
5x, 


^ (x) 
5x7 


5f 2 (x) 
6x„ 


afj (x) 
5x~ 


5f 2 (x) 
5x_ 


8f n (x) 
Sx_ 


( A— 7 ) 


(X) 


<V 


Next (F(X q )} is evaluated by using Eq. (A-6). With these values 
°f (F(X o )} and [J(X Q )] the next approximation X^ is obtained as. 


tX l> ■ <V„xl - [W < X o>C (F < X o»n + l] ,fl - 8) 


With this new value of {X^, (F( X ^ ) > and [ J ( X-^ ) ] are calculated and 

then { X 2 > is calculated as: 


<V ' < X l>nxl - «1»C (F < X l»nxll 


( A— 9 ) 


This iterative procedure is continued till the maximum difference 
in (X) values between the successive iterations is within specified 
tolerence of e i.e. 


Max j (X i > - tX i -_ 1 } | < e ( A— 10 ) 

When Eq. (A-10) is satisfied the iterations are stopped and the 
solution of Eq. (A-5) is {X.}. 
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A. 3 Example Problem 

In the illustrative example (section 3.3.2), a circular membrane 
with radius a and thickness h, and with a uniform loading over the 
region 0 < — < 0.5 was analyzed. By using the four region idealiza- 
tion and a finite difference technique, five nonlinear equations were 
obtained. These five nonlinear equations (Eqs. (3.26), (3.27), (3.28), 
(3.29) and (3.24)) were expressed in terms of normalized stresses as: 


- 3 a.(0) + 4 0.(1) - o r (2) = 0 ( A- 1 1 ) 

- 0.5 3 .( 0 ) o. 2 (l) - 2 o. 3 (l) + 2.5 0.(2) o. 2 (l) = - ^ (A-12) 

5(1) 5 2 (2) - 8 0.3(2) + 7 o (3) o 2 (2) = - Jr (A-13) 

r r r r r 3z 

4.5 0.(4) o r 2 (3) + 1.5 0.(2) o. 2 (3) - 6 o. 3 (3) = - (A-14) 

0,(4) (6.69) - 8 o.(3) + 2 0.(2) = 0 (A-15) 


where 



Using Newton-Raphson method [J (o^,)] was obtained as: 
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-3 4 -1 

- o p (0) 0,(1) 

-0.5 o. 2 ( 1) - 6 o. 2 ( 1) 2.5 o 2 (l) 

+5 a (2) o (1) 
r r 


0 0 
0 0 


0 o r 2 (2) 2 0.(1) o.<2) 7 c. 2 (2) 0 

- 24 a r 2 (2) 

+ 14 o r (3) 0.(2) 

0 0 1.5 o r 2 (3) 9 o f (4) o p (3) 4.5 a r 2 (3)j 

+ 3 0.(2) 0.(3) 

- 18 o 2 (3) 

0 0 2 - 8 6.69 _j 

(A -16) 


An initial approximation for the normalized stress 0^(0), 0^.(1), 

0^.(2), 0^(3), and 0^(4} was assumed as unity. With this assumption, 
the function (F(o r )} and [J( 0 r )] were evlauated. They are 


' 

0 


- 3 

4 

-1 

0 

0 

1/128 

_ 

-0.5 

-2 

2.5 

0 

0 

1/32 

It 

l — l 

J 

1 t) 

\ I 

0 

1 

-8 

7 

0 

1/216 

r 

0 

0 

1.5 

-6 

4.5 

0.69 


0 

0 

2 

-8 

6.69 


( A— 1 7 ) 


The next approximation for (o } was obtained by using Eq. (A-8) as: 

f , /■ 


o r (0) 


/■ - 
1 


I 

-3 

4 

-1 

0 

0 


/ -\ 
0 

* r d) 


1 

- 

-0.5 

-2 

2.5 

0 

0 


1/128 

5 r< 2) 

■ { 

1 


0 

1 

-8 

7 

0 

( 

1/32 

W««‘ 


1 


0 

0 

1.5 

-6 

4.5 


1/216 

4* 


1 

) 


0 

0 

2 

-8 

6.69 


0.69 

(A-18! 
< >/ 
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Solution of Eq. (A-18) is 


N 



a r (0) 


0.02039 

0.(1) 


0.01941 

a (2) 

r 

■ 

0.01648 

0.(3) 


0.01160 

o (4) 


0.0089 

k r , 




The above procedure was repeated and the 

to (a } was obtained as: 
r 


( A— 19 ) 


next approximation 


/ >j 


r- 

o p (0) 


27.5628 

o r (l) 


27.7630 

o r (2) 


28.3637 

o r (3) 


27.3757 

0.(4) 


24.2568 

k J 


J 


This iterative procedure was continued till the maximum difference 

in (a } values between the successive iterations was less than IE-6 or 
r 

max | A o r ' j < IE-6 . The solution converged after 18 iterations, 

and the corresponding normalized stress values were: (0) = 0.3156, 

a (1) = 0.3065, 5 (2) = 0.2726, 5 (3) = 0.2233, and S (4) = 0.1851. 
p r r r 
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APPENDIX B 

LARGE DEFLECTION SOLUTION OF CLAMPED CIRCULAR PLATE LOADED AT 
THE CENTER USING THE ENERGY METHOD 

Consider clamped circular plate of thickness h and radius a, 
subjected to a center point load P. The classical large deflection 
solution based on the energy method assumes that the deformation shape 
of a clamped circualr plate under center point load has the same 
equation as in the case of small deflections, therefore 

w(r) = w q [1 - ( 4 ) + 2 ( 4 ) (!) ] (B-l) 


where w Q is the central displacement of the plate. 

The corresponding strain energy of bending can be written as [27] 


D 2n a ,2 0 , jo 2v . ,2 , 

u b = -|a / / [(^) 2 + 4(^} 2 + -|5^i L^rdrde 

o o dr r dr 


( B— 2 ) 


8 7i D w 
eq o 


( B-3) 


where D eq is a flexural modulus and is given by 


D 


E h~ 
eq 


eq 12 (l-v\) 
eq 


( B— 4 ) 
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where E gq and v gq are the equivalent Young's modulus and Poisson's 
ratio for the quasi-i sotrpic plate and were obtained by using the 
technique described in Chap. 2. 

For the radial displacement u, classical solution assumes the 
expression 


u = r (a - r) (C^ + C 2 r) 


( B— 5 ) 


Equation (B-5) satisfies the boundary conditions that u must vanish 
at the center and at the edge of the plate. 

The strain energy due to stretching of the middle plane of the 
plate is given by [27] 

u m - 77H7 * <4 * 4 * 2 % V E e> r dr (B - 6) 

(1_v eq ) 0 

where e r and e Q are radial and tangential strains and are given by 
Eqs (3.3) and (3.4) from Chap. 3 as: 




r dw ^ i 


( B— 7 ) 


E 


9 



( B— 8 ) 


From expressions (B-l) and (B-5) for the displacements, the strain 
components e r and e Q were calculated by using Eq. (B-7) and (B-8). 
Substituting these strain components in Eq. (B-6) strain energy due to 
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stretching U m was obtained as 


U = 
m 


71 ^ ha o 5 ? . 

rS — 5— [0.250 C7 a^ + 0.1167 C, a* 

» - v eq> 


w 2 


+ 0.300 C. C 0 a 3 + 0.06332 C. 

12 la 


( B— 9 ) 


w 


+ 0.09851 C 0 + 0.18240 -Z-] 
Z 0 _4 


The constants Cj and C 2 were determined from the condition that the 
total energy of the plate for a position of equilibrium is minimum. 
Hence 


8U 


dU 


aq = 0 and aq ~ 0 


( B— 10 ) 


Substituting Eq. (B-9) for U m , following two linear equations for 
Cj and C 2 were obtained 


w 

0.50 C. a 2 + 0.300 C 9 a 3 - 0.06332 — 

1 c. 3 


( B— 1 1 ) 


0.2334 C, a 4 + 0.300 C. a 3 = - 0.09851 w 2 (B-12) 

2 1 0 


Equations (8-11) and (B-12) were solved simultaneously and and 
C 2 were obtained as: 


2 2 
w w 

C. = 0.5531 = -1.133 

1 a 3 2 a 4 


(B-13) 



190 


Substituting values of Cj and C 2 from Eq. (B-13) into Eq. (B-9), U m 
was obtained as 

4 

( B— 14 ) 


w_ 


11 h 

U_ = - q - 7 ■ [0.14409 -^] 


m 


1 1 - V 1 


Total potential energy n was written as the sum of the bending energy 
and stretching energy minus the work done, therefore 


8 it D w 


w H it E h 
eq 


■ 4 3 —°- + 0.114409 JL r 

a a (1 - v „ ) 

eq 


- P w 


( B— 1 5 ) 


Minimizing iip with respect to w Q , following relationship between load 
P and w Q was obtained 


w w it E h 

16 it D -4 + 0.57636 -2 ™ — 

eq ^ ^ (1 - v eq ) 


= p 


( B- 16 ) 


E eq 


Substituting D = — 

eq 12 (1 - V 
• G Q 

can be written as: M 


Y~ . and rearranging the terms, Eq. (B-16) 


w w 3 

(-ft) + 0.433 (—ft) = 0.2157 


pa 


E eq h 


( B— 1 7 ) 


Equation ( B- 1 7 ) is the classical large deflection solution for 


clamped circular plate under central point load P. 
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APPENDIX C 

DERIVATION FOR THE RADIUS OF INFLECTION 

The classical large deformation theory [27], assumes the deforma- 
tion shape of a clamped circular plate under central point load as: 


w(r) 


w 


[1 - (^) + 2 (ly) An 

a 


a 




(C-l) 


Furthermore the classical theory assumes that the functional form of the 

deformation shape remains unchanged for various values of central 

deflection, w Q . Thus the radius of inflection (the radius at which 

w 

curvature is zero) is independent of (— . The radius of inflection 
can be obtained as follows: 

By differentiating Eq. (C-l) twice, curvature can be written as 

( C— 2 ) 

(^) - 0 . 
dr 

( C— 3 ) 


.2 4 w 

d w o 


[1 + An [L)] 


dr a 

The radius of inflection is the radius at which 
Hence, 

[1 + {[)] = 0 


yields the radius of point of inflection as 

(L) = e" 1 = 0.3678 

a 


(C-4) 
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